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ABSTRACT 


Recent  interest  in  pollution  control  and  the  proximity  of  Monterey 
Bay  to  the  Naval  Postgraduate  School  prompted  an  investigation  of  the 
circulation  in  the  bay.   The  first  phase  of  the  study  consists  of 
solving  the  simple  cavity  flow  problem.   A  vorticity-stream  function 
relationship  is  solved  using  an  explicit,  time  dependent,  finite  dif- 
ference scheme.   Solutions  for  selected  Reynolds'  numbers  and  length 
to  width  ratios  of  the  cavity  are  obtained.   Values  are  chosen  to  give 
an  indication  of  the  flow  patterns  occurring  over  a  wide  range  of 
these  parameters. 

Equations  for  a  refined  model  are  derived  to  include  the  effects  of 
the  bottom  topography,  frictional  forces  and  the  Coriolis  force.   A 
numerical  procedure  similar  to  the  one  applied  to  the  simple  cavity 
flow  problem  is  used  on  the  refined  equations.   The  topography  of 
Monterey  Bay  is  used  in  this  study.   Results  indicate  that  closed 
circulation  patterns  are  more  probable  in  Monterey  Bay  than  in  other 
embayments  due  to  the  presence  of  the  submarine  canyon. 
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I .   INTRODUCTION 

A.   BACKGROUND 

Monterey  Bay,  located  on  the  West  Coast  of  the  United  States 
approximately  one  hundred  and  twenty  miles  south  of  San  Francisco, 
presents  a  good  opportunity  to  study  currents  within  an  embayment 
(Figure  1.1).   This  bay  is  well  suited  for  circulation  studies  because 
of  the  following  geographic  and  oceanographic  features.   The  bay  is  a 
large  one  and,  to  a  close  approximation,  it  has  symmetrical  boundaries. 
It  is  semi-elliptical  in  shape  with  a  major  axis  of  twenty  miles  and 
a  width  of  about  eight  miles  (Figure  1.2).   Bisecting  the  bay  is  a 
large  submarine  canyon.   This  canyon  originates  in  the  bay  offshore  of 
Moss  Landing.   At  the  mouth  of  the  bay,  eight  miles  from  Moss  Landing, 
the  canyon  attains  a  depth  of  six-thousand  feet.   This  topographic 
feature  is  as  large  as  the  Grand  Canyon  in  vertical  relief.   The  effect 
of  the  Monterey  Submarine  Canyon  on  the  circulation  within  the  bay  is 
one  of  the  problems  which  is  investigated  in  this  study.   A  desirable 
oceanographic  feature  is  the  presence  of  currents  which  flow  past  the 
bay.   The  Davidson  Current  flows  north  along  the  West  Coast  of  North 
America  from  November  to  February.   Flowing  southward  during  the 
spring  and  summer  is  the  California  Current.   It  is  assumed  in  this 
study  that  oceanic  currents,  such  as  these,  are  the  driving  force  for 
the  circulation  within  an  embayment.   The  currents  drive  the  circulation 
by  advecting  momentum  into  the  area  in  question.   The  effects  of  tidal 
forces  and  local  winds  are  neglected.   It  is  not  meant  to  be  implied 
that  these  forces  are  unimportant,  but  that  only  a  single  component  of 
the  total  system  is  being  considered  in  this  analysis. 
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Monterey  Bay  is  conveniently  located  for  conducting  circulation 
studies  for  several  reasons.   First,  the  Naval  Postgraduate  School  is 
located  on  the  shores  of  the  bay  in  the  city  of  Monterey.   Both  the 
Postgraduate  School  and  the  United  States  Navy  are  interested  in  circu- 
lation studies  of  this  kind.   The  facilities  at  this  institution  are 
well  adapted  to  carry  out  such  studies.   The  Department  of  Oceanography 
at  the  Naval  Postgraduate  School  is  interested  in  circulation  studies 
and  has  access  to  oceanographic  research  vessels.   Also,  located  at  the 
School  is  a  large  computer  center  with  a  well-trained  staff  and  an 
IBM  360/67  Computer  along  with  all  the  necessary  peripheral  equipment. 
Secondly,  there  are  many  communities  located  on  the  shores  of  Monterey 
Bay  which  are  becoming  increasingly  aware  of  pollution  problems.   An 
important  question  that  needs  to  be  answered  before  pollution  can  be 
controlled  is:   what  is  the  effect  of  currents  on  the  discharged 
effluent?   To  answer  this  question  the  circulation  regime  of  the  area 
involved  must  be  determined.   It  is  hoped  that  the  results  of  this  study 
will  provide  some  answers  to  the  above  question. 

Prior  to  this  study,  little  had  been  done  to  ascertain  the  circu- 
lation patterns  within  Monterey  Bay.   This  thesis  is  part  of  a  larger 
research  effort  being  conducted  by  the  Naval  Postgraduate  School  aimed 
at  determining  the  forces  which  dominate  in  controlling  current  patterns 
within  an  embayment  and  the  types  of  patterns  which  are  generated.   Pre- 
sently, there  are  three  phases  to  the  project.   Two  phases  concern 
numerical  modeling.   One  numerical  model,  the  subject  of  this  thesis, 
assumes  oceanic  currents  drive  the  circulation  within  a  bay.   A  second 
model  simulates  the  tide  and  wind  induced  circulation.   The  third  phase 
involves  field  studies  conducted  in  order  to  assess  the  models  which 
have  been  developed. 
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Figure  1.1  Location  of  Monterey  Bay 
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Figure   1.2     Monterey  Bay 
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B.   NUMERICAL  MODELING 

Recently,  with  new  developments  in  partial  differential  equation 
theory  and  finite  difference  methods,  the  field  of  fluid  mechanics  has 
progressed  tremendously.   The  advent  of  large-scale  computers  has 
spurred  these  developments  and  made  it  possible  to  solve  partial  dif- 
ferential equations  numerically  which  have  proven,  so  far,  to  be 
impossible  to  solve  analytically.   Due  to  the  great  speed  of  these 
computers  it  has  become  possible  to  solve  large  systems  of  simultaneous 
equations  rapidly  enough  to  put  numerical  prediction  models  in  fields 
such  as  Meteorology  and  Oceanography  on  a  real  time  basis. 

Most  circulation  problems,  basically,  involve  the  solution  of  a 
large  number  of  simultaneous  equations  with  specified  boundary  con- 
ditions.  Practically  all  studies  of  this  type  are  concerned  with 
either  small-scale  flows  with  simple  geometries  such  as:   channel  flow, 
flow  around  a  cylinder,  and  simple  cavity  flow  or  with  large-scale  oceanic 
circulation  studies.   Little  has  been  done  to  try  to  apply  these  equations 
and  numerical  techniques  to  medium-scale  flow  such  as  the  circulation 
within  an  embayment.   There  are  several  possible  reasons  why  more  medium- 
scale  flow  problems  have  not  been  investigated.   One  reason  is  the  high 
cost  of  studies  of  this  type.   Computer  time  is  very  expensive  and  not 
always  readily  available.   Medium-scale  problems  are  in  general  highly 
variable  and  each  case  must  be  dealt  with  individually.   Another  reason 
is  that,  before  the  increased  awareness  of  pollution  problems,  there  was 
little  interest  in  the  regional  circulation  beyond  that  caused  by  waves 
and  tides.   Lastly,  the  application  of  the  equations  of  motion  to  large 
or  medium-scale  flow  problems  requires  assumptions  about  the  forces  and 
terms  involved  which  have  not  yet  been  universally  accepted  as  being  valid, 
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There  are  many  numerical  procedures  that  can  be  used  to  try  to  solve 
a  system  of  simultaneous  partial  differential  equations.   An  interesting 
result  of  applying  different  numerical  techniques  to  the  same  problem  is 
the  possibility  of  obtaining  different  solutions.   One  of  the  problems  of 
numerical  modeling  is  determining  whether  the  correct  solution  has  been 
reached,  if  only  one  exists.   There  has  been  a  lot  of  research  into  the 
problem  of  the  solution  of  partial  differential  equations  by  Ames  [Ref.  2], 
Considering  the  Navier-Stokes  Equations,  he  proposes  that  a  unique 
solution  exists  only  below  a  certain  Reynolds'  number  and  that  above  this 
Reynolds'  number  several  solutions  exist;  whereas,  above  a  larger'  critical 
number  no  solutions  exist.   The  critical  Reynolds'  number  marks  the 
transition  to  turbulent  flow.   To  illustrate  the  existence  of  multiple 
solutions  to  a  partial  differential  equation,  he  gives  an  example  of  a 
quasi-linear,  elliptical  partial  differential  equation  for  which  there 
is  no  unique  solution. 

An  extensive  study  of  various  numerical  procedures  used  in  the 
solution  of  compressible  and  incompressible  laminar  separated  flow  for 
simple  geometries  was  made  by  Roach  [Ref.  7].   His  dissertation  is  quite 
thorough  in  its  consideration  of  the  problems  encountered  and  the 
techniques  used  to  resolve  them.   Some  of  Roach's  techniques  are 
extended  in  this  study  to  the  problem  of  motion  within  an  embayment  in 
the  simple  cavity  flow  section  of  this  thesis. 

C.   SCOPE  OF  PRESENT  PAPER 

This  thesis  attempts  to  take  the  equations  of  motion  and  apply  them 
to  the  problem  of  the  circulation  in  an  embayment  and  solve  them 
numerically.   A  progression  is  made  from  the  simple  cavity  flow  problem 
to  the  refined  model.   The  simple  cavity  flow  problem  considers:   local 
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rate  of  change  of  vorticity  with  time,  advective  terms,  and  lateral 
shear  stresses.   In  addition,  the  refined  model  considers:   planetary 
vorticity  tendency,  topographic  vorticity  tendency  and  bottom  friction. 
Procedures  for  numerically  solving  these  equations  are  examined  while 
economy  of  computer  time  and  the  validity  of  the  solution  are  kept  in 
mind.   An  attempt  is  made  to  generalize  the  computer  program  so  that  it 
can  be  used  on  any  embayment.   The  generality  is  obtained  by  allowing  the 
input  of  an  arbitrary  bottom  topography. 

The  model  is  tested  on  Monterey  Bay  in  order  to  discover,  not  only 
the  circulation  patterns  resulting  from  the  interaction  of  the  forces 
involved,  but  also  the  relative  importance  of  the  forces.   To  determine 
the  importance  of  the  forces  involved,  various  computer  runs  are  made 
with  different  combinations  of  the  forces  acting,  and  the  effect  on  the 
circulation  patterns  is  examined.   It  is  hoped  that  the  results  of  these 
investigations  will  shed  some  light  on  the  processes  which  control  circu- 
lation patterns  and  help  provide  answers  to  some  of  the  problems  which 
now  exist. 
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II.      NUMERICAL  METHODS 

A.  INTRODUCTION 

This  section  consists  of  a  brief  discussion  of  some  of  the  aspects 
of  numerical  modeling  which  are  used  in  this  thesis.   There  are  many 
excellent  texts  which  contain  thorough  discussions  of  this  material; 
therefore,  the  material  is  not  covered  in  depth.   All  of  the  repre- 
sentations derived  in  the  following  section,  II-B,  will  be  used  in  the 
development  of  the  numerical  equations  solved  in  later  sections.   The 
following  paragraphs  introduce  these  concepts  and  develop  the  necessary 

relationships . 

♦ 

B.  NUMERICAL  RELATIONSHIPS 

To  represent  a  partial  differential  equation  in  numerical  form  the 
starting  point  is  usually  Taylor's  Theorem.   From  this  theorem,  the 
necessary  expressions  which  represent  derivatives  in  their  various  forms 
can  be  derived.   Numerical  relationships  for  terms  such  as  the  Laplacian 
and  the  Jacobian  can  then  be  developed  from  these  expressions.   Basically, 
these  expressions  determine  the  value  of  a  derivative  at  a  given  point 
in  terms  of  the  values  at  surrounding  points  and  the  grid  spacing. 
1.   Taylor's  Theorem 

The  finite  difference  representation  for  the  derivatives  of  a 
function  at  a  given  point  in  terms  of  the  values  of  the  function  at 
surrounding  points  can  be  derived  from  Taylor's  Theorem. 
Taylor's  Theorem  in  Two  Dimensional  Space: 

If  u  and  its  derivatives  of  order  _<  p+1  are  single  valued, 
finite  and  continuous  at  every  point  on  the  line  segment 
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from  (x,y)  to  (x+Ax,y+Ay)  then  there  is  a  point 
(a,b)  on  that  line  segment  such  that 

P         3N+M 
u(x+Ax,y+Ay)  =  u(x,y)  +    E    ^77 '"f^^l    (Ax)N  (Ay)M  +  R    (2.1) 

(N+M)=l   '  *  3x3  y 


where  the  remainder,  R,  is  given  by 

9N+M 
t.  !  v  u(a,b)     -      .N    ,.    .M 

R  =  rar  N+M  +,~J^  (  }   ( y) 

N+M=p+1  3x    dy 

It  can  be  seen  that  if  Ax  and  Ay  <  1  as  N  and  M  ■*■  «  then  R  -*■   0 ; 
and  in  the  limit,  the  expansion  equals  the  function.   If  y  is  held 
constant  (i.e.,  Ay  =  0)  there  is  a  contribution  to  the  summation  terms 
only  when  M  =  0,  therefore 


P      N 
u(x+Ax,y)  =  u(x,y)  +  £  |t  ^~  (Ax)N  +  0  (xP+1)  (2.2) 

N=l   *  3x 


If  x  is  held   constant,    (Ax  =  0),    similarly 


P  M 

u(x,y+Ay)    =  u(x,y)   +     I     jfr  ^  (Ay)M  +  0    (yP+1)  (2.3) 

M=l  3y 


Using  these  equations,  expressions  for  derivatives  at  a  point 
are  developed  in  Section  II-B-2. 

2.   Representation  of  Derivatives 

The  following  grid  illustrates,  diagrammatical ly ,  the 
notation  that  is  used  to  fix  the  grid  position  of  a  given  point. 
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x-Ax,y+Ay 


x,y+Ay 


x+Ax , y+Ay 


Ax 


Ui-l,j+l 


x-Ax,y 


Ay 


l.j+l 


x,y 


i+l,j+l 


x,y+Ay 


u. 


1-l.J 

x-Ax,y~Ay 


u 


x,y-Ay 


um,j 


x+Ax,y-Ay 


i-l.j-l 


Ui,j"l 


1+l.J-l 


Using  Taylor's  Theorem  with  y  held  constant  Equation  2.2  in 
the  new  notation  becomes 

3u, 


_LLL 


a2 

0  u.   . 


u;x1  ,  ■  u.  .  +  — ^-   Ax  +  , 
l+l,  j    i,j     3x       2   3x2 


1    1,1  -  .2  .  .  ,.  .3 


(Ax)   +  0  (Ax) 


If  -  Ax  is  substituted  for  Ax  the  equation  becomes 

9u. 


32u 


u.  -  .  =  u.  .  -  — ±-^  Ax  +  ±  i^l  (Ax)Z  +  0  (Ax) J 

l-l, j    i,j     3x        2    3x2 

Linear  combinations  of  2.4  and  2.5  with  the  higher  order  terms 
neglected  yields 

Representation  Mag,  of  Error        Diff.  Scheme 


a    u .-  .  -  u.  . 

3x         Ax 


u.  .  -  u     . 


ax 


Ax 


0  (Ax) 


0  (Ax) 


Forward 


Backward 


(2.4) 


(2.5) 


»      u,  n   .  -  u .  n 

3u  =   i+l, ,1     i-l,,1 

9x        2Ax 

.2    u...  .  +  u.  1  . 
3  u  _   i+l,, 1     i-l,, 1 


-  2u, 


3x 


(Ax) 


0  (Ax) 


i*i     0  (Ax): 


Central 
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Similarly,  derivatives  with  respect  to  Y  can  be  found. 

Representation  Mag,  of  Error  Diff.  Scheme 


a     u-  -.1.1  "  u-  • 
3y         Ay 


r.  U.    .   "   U. 


Ay 


0  (Ay)  Forward 

0  (Ay)  Backward 


3y        2Ay 


0  (Ay)  Central 


.2  u.  ,   +  u.  .  _  -  2u.  . 

A-"  =  -Uii ^J"1 i-J.       n  M,^2 

2  2 

3y  (Ay)Z 


0  (Ay)' 


A  representation  for  the  Laplacian,  V  u,  can  be  developed 
from  the  above  expressions. 


2     2 

„2    3  u  .  3  u 
V  u  =  — 7?   + 


2     2 

3x    3y 


o          u.,n     ,   +  u.    .     .    -   2u.    .        u.     . ,.   +  u.     .    -    -   2u. 
y2u  =     1+1>J 1~1>.1 JiX  +     1?Jtl 1>.l-1 LA  (2,6) 

(Ax)2  (Ay)2 

3.   Representation  of  the  Arakawa  Jacobian 

This  study  utilizes  Arakawa' s  technique  for  representing  the 
Jacobian  which  is  an  average  of  three  different  expressions  for  the 
Jacobian  [Ref.  3].   The  method  Arakawa  developed  is  used  for  two 
reasons.   First,  it  has  the  desirable  property  of  conserving  certain 
quantities  two  of  which  are  the  vorticity  and  the  vorticity  squared. 
A  further  discussion  of  the  conservative  property  is  contained  in 
Section  II-D.   Secondly,  the  representation  is  a  very  stable  one  which 
helps  to  minimize  convergence  problems  in  solving  the  numerical  equations, 
The  price  that  is  paid  for  this  stability  is  the  high  degree  of  smoothing 
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inherent  in  this  representation  which  may  eliminate  some  small-scale 
features.   Although  these  features  may  be  real,  it  is  not  felt  that 
their  elimination  is  very  important  because  this  study  is  only  inter- 
ested in  the  general  flow  patterns.   Examination  of  small-scale  features 
in  the  flow  patterns  with  the  simplifying  assumptions  which  are  made  in 
this  study  would  not  be  realistic. 

In  order  to  develop  the  desired  relationship  for  the  Arakawa 
Jacobian,  the  following  three  expressions  for  the  Jacobian  are 
considered. 

...    _.    9A  3B    dk   3B    TA  .        . 

J(A,B)  =  — —  —  =  J  (2.7) 

3x  3y    dy  dx 

J<A-B>  ■  h  <A  !f>  -  k  (A  £ }  =  jB  (2-8) 

J<A>B>  -  i  <B  $  -  h  <B  § >  - jC  <J-9) 

The  Arakawa  Jacobian  is  given  by 

j  =  §  (/  +  jb  +  jc) 

Using  the  same  notation  as  in  the  previous  section,  the  expressions 
become 

jA  =  — L   [(A.,-  .-A.  .  .)(B.  ..--B.    .)  -  (A.  ...-A.  .  tXB,,-  .-B.  .  .) 
4AxAy    i+1,3   i-l, j   i,j+l  i,j~l      i,3+l   i,j~l   i+l»J   i-l»3 

JB  =  rr-h      [A_  .  (B...  .,,-B'    .  -)  -  A.  ,  .  (B.  .  .,,-B.  -  .  .) 
4AxAy   i+l,3   i+l,j+l  i+l,j-l     i-l,J    i-l.j+l  i-l>3"l 

1,3+1   1+1,3+1   i-l,3+l     1,3-1   1+1,3-1   1-1,3-1 

JC  ■  ;  \      [B.  ..-(A..-  ,.,-A..-  .  -)  -  B.  .  -  (A.  n  ,..,-A.  ,  ,  J 
4AxAy  L  i,3+l  1+1,3+1  1+1,3-1     i,j-l   1-1,3+1  i-l,j-l 

-  B...  .(A.,-  ,.,-A,.,  .  .)  +  B.  ,  .(A.  .  ,.,-A.  .  .  ,)] 
1+1,3   1+1,3+1   1+1,3-1     i-l,3   i-l,3+l   1-1,3-1 
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(2.10) 


Collecting  terms  and  factoring  yields  the  desired  form. 


J(A,B)  =   l  [B...  .  (A.  .  ,+A.^  .   -A.  ..--A..-  ..-)  + 

12AxAy   i+l,3    i,j-l   1+l.j-l   i,j+l   1+1, j+1' 

B4   «X1   (A--l-1   -+A--i.1   -U-l-A-   1   -_A-   1   iA.0     + 

i,j+l   i+l,  J   i+l,  3+1   i-l,J   1-1,3+1 

B.  .  .  (A.  .,,+A.    ,,,-A.  .   -A.    .   )  + 
i-l, J    1,3+1   1-1,3+1   1,3-1   1-1,3-1 

B.  .  .  (A.    ,+A.    .  --A..-  ,-A...  .  .)  + 
1,3-1   i-l,3   1-1,3-1   i+l, J   1+1,3-1 

1+1,3+1   i+l,3   1,3+1 

B.  .  .,-  (A.  ,,,-A.  .  .)  + 
1-1,3+1    1,3+1   1-1,3 

B.    .    (A.    .-A.  .  J  + 
1-1,3-1    i-l,3   1,3-1 

B...  .  -  (A.  .  ,-A..,  .)] 

i+l, 3-1   1,3-1   1+1,3 


C.   RELAXATION  METHODS 

One  general  approach  to  solving  a  set  of  simultaneous  equations  is 
the  relaxation  method.   This  method  attacks  the  problem  by  determining 
the  difference  between  an  estimated  solution  to  the  equations  and  the 
actual  solution.   Each  successive  solution  to  the  equations  differs 
from  the  "true"  solution  by  a  determined  amount.   The  aim  of  this  method 
is  to  manipulate  the  equations  so  that  the  differences  from  the  "true" 
solution  approach  zero.   The  manner  in  which  this  is  achieved  is 
explained  thoroughly  in  Ref.  1  and  discussed  briefly  in  the  following 
three  sections,  (II-C-1,  II-C-2,  II-C-3) .   This  discussion  begins  with 
the  definition  of  a  residual. 

1.   Residual 

One  can  write  a  set  of  simultaneous  equations  such  that  all  the 
terms  are  on  one  side  of  the  equality  sign,  for  example: 
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X  +  Y  -  3  =  0 

Residuals  are  defined  as  quantities  which  represent  the  values  of  the 
equations  when  values  are  assigned  to  the  variables,  X  and  Y.   The 
solution  to  the  equations  is  reached  when  the  residuals  equal  zero. 

X  +  Y  -  3  =  R  (2.11) 

x 

2X  +  3Y  +  2  =  R  (2.12) 

y 

The  residuals,  R  and  R  ,  are  a  measure  of  the  closeness  of  X  and  Y  to 
x      y 

the  exact  solution.   The  problem  of  solving  simultaneous  equations 

then  becomes  one  of  finding  values  for  X  and  Y  which  make  the  residuals 

approach  zero.   One  method  for  doing  this  is  called  relaxation. 

2.   Relaxation 

The  relaxation  method  involves  taking  the  largest  residual  and 

computing  a  value  for  its  associated  variable  which  makes  the  residual 

zero.   Using  the  above  two  equations  as  an  example,  the  variable  X  is 

arbitrarily  associated  with  the  residual  in  Equation  2.11.   Likewise,  the 

variable  Y  is  associated  with  the  residual  in  Equation  2.12. 

Assuming  that  R  is  the  larger  residual,  then  to  make  R  =  0 

only  the  value  of  X  would  be  adjusted.   Now  that  R  =  0,  R  is  the 

j  x       y 

largest  residual,  and  the  value  of  Y  in  Equation  2.12  would  be  adjusted 

to  make  R  =0.   This  process  is  continued  until  the  values  of  R  and  R 
y  x      y 

are  within  acceptable  limits. 

It  is  not  economical  in  terms  of  computer  time  to  determine 
which  equation  has  the  largest  residual  and  relax  only  it.   Making  many 
sweeps  through  all  of  the  equations  and  relaxing  them  simultaneously  is 
a  more  efficient  method.   The  reason  for  this  is  that  logical  statements 
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take  much  more  computer  time  to  execute  than  arithmetic  statements. 
The  process  which  relaxes  all  of  the  equations  simultaneously  until 
all  of  the  residuals  are  smaller  than  a  predetermined  value  is  called 
simultaneous  relaxation.   Modifying  this  procedure,  it  is  found  that 
even  more  computer  time  is  saved  if  simultaneous  overrelaxation  is 
used. 

3.   Overrelaxation 

Overrelaxation  is  a  process  whereby  the  residuals  are  not 
reduced  to  zero  but  are  overcompensated  for  in  order  to  speed  the  rate 
of  convergence.   This  is  the  usual  process  used  because  it  minimizes 
the  computer  time  needed  to  solve  the  equations.   This  method  is 
approached  by  writing  the  numerical  equations  in  such  a  way  that  the 
residuals  are  multiplied  by  a  convergence  factor  which  is  called  the 
optimum  overrelaxation  factor,  Ropt.   A  critical  part  of  the  problem 
is  determining  the  factor  which  will  yield  the  fastest  convergence  rate. 

A .   Optimum  Overrelaxation  Factor,  Ropt 

a.  Introduction 

A  factor,  Ropt,  can  be  determined  which  maximizes  the 
convergence  rate  of  a  set  of  simultaneous  equations.   The  convergence 
factor  must  be  between  1.0  and  2.0  and  its  choice  may  be  critical, 
because  small  differences  in  the  factor  may  make  large  differences  in 
the  convergence  rate.   There  are  two  ways  of  determining  Ropt.   One 
way  is  through  theoretical  considerations  and  the  other  is  to  determine 
it  empirically.   These  methods  are  discussed  in  the  following  two 
sections. 

b.  Theoretical  Ropt 

Application  of  theoretical  considerations  must  be  restricted 
to  simple  partial  differential  equations  and  simple  boundaries  because 
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of  the  complexity  of  the  problem  [Ref.  7].   Frankel  applied  the  concept 
of  overrelaxation  to  Leibmann's  method  for  solving  the  following 
equation 

Z  =  V2^  (2.13) 

Leibmann's  method  for  solving  the  above  equation  involves  deriving  its 
numerical  representation  and  solving  for  H*.   This  method  produces 


yk+1  m     k     _R££t   [VFk     +  M       +  B2   k     +  ^k+1 
i»j    i,J   2(1+B  )    1+1»1    1-1»J       '  i»J+l    i,J-l 

9   k  ?    k 

-  (Ax)   Z   ,  -  2  (1  +  B  )  V   .] 


(2.14) 


The  bracketed  term  represents  the  residual,  and  the  solution  is  being 
reached  as  it  approaches  zero. 

Frankel  determined  that  for  a  rectangular  grid  with 
Dirichlet  boundary  conditions  the  optimum  overrelaxation  factor  is 
given  by 


where 


and 


1  _  A-a 
Ropt  =  2  [- v-±-±]  (2.15) 


cos  (£)  +  B   cos  (J) 

a  =  [ 3 —] 

1  +  B 


Ax 


B  - 

Ay 


N  =  Number  of  rows  (Ax) 
M  =  Number  of  columns  (Ay) 
It  is  found  that  for  a  41  x  80  grid: 
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B  =  1    Ropt  -  1.8574 
B  =  2    Ropt  =  1.8916 
B  =  3    Ropt  =  1.9064 
c.   Empirically  Determined  Value 

It  can  be  seen  from  the  last  section  that  Ropt,  theoretically, 
depends  upon  the  grid  size  ratio  as  represented  by  B  and  the  dimensions 
of  the  array  being  relaxed.   A  study  was  conducted  to  experimentally 
determine  Ropt  for  the  following  grid  which  is  the  one  used  in  the 
simple  cavity  flow  problem. 


41    GRID    POINTS 


80 


10 


21 


21 


10 


Figure  2.1  Numerical  grid 

Leibmann's  equation  for  solving  for  the  stream  function  is  the  one 
relaxed.   This  is  Equation  2.14  with  Ropt  =  1. 

Evaluating  the  equation  involves  initializing  the  grid  and 
relaxing  the  equations  until  every  value  on  the  grid  differs  from  the 
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previous  value  by  less  than  1.0  x  10   .   This  is  approximately  the 
limit  of  accuracy  on  the  IBM  360/67  Computer  without  resorting  to 
double-precision  accuracy  which  would  significantly  increase  compu- 
tational time.   To  further  optimize  use  of  computer  time  the  convergence 
is  checked  only  once  every  ten  passes  through  the  grid.   Each  pass  con- 
sists of  a  sweep  through  the  grid  from  left  to  right,  top  to  bottom 
followed  by  a  sweep  from  right  to  left,  bottom  to  top.   A  limit  of 
two-hundred  passes  through  the  grid  is  placed  on  each  test. 

Results  of  this  study  are  summarized  in  Table  I   and  pre- 
sented in  graphical  form  in  Figures  2.2,  2.3,  and  2.4.   Table  I   gives 
the  number  of  passes  through  the  grid  required  for  every  value  on  the 
grid  to  converge.   The  criterion  for  convergence  is  that  the  difference 
between  two  successive  passes  is  less  than  1.0  x  10   .   Results  are 
given  for  four  different  Reynolds'  numbers  and  for  aspect  ratios  of 
1,  2,  and  3  represented  by  Ax  =  .05,  .10,  and  .15  respectively.   In  all 
cases  Ay  is  held  constant  at  .05.   Figures  2.2,  2.3,  and  2.4  represent 
the  same  data  in  graphical  form.   This  illustrates  the  effect  of  Ropt 
on  the  convergence  of  the  stream  function  equation  in  a  more  explicit 
manner . 

One  conclusion  which  can  immediately  be  drawn  from  these 
empirical  results  is  that  Ropt  is  also  a  function  of  the  Reynolds' 
number.   The  most  significant  result  of  this  study  is  that  overestimating 
the  value  to  be  used  for  Ropt  may  significantly  increase  the  convergence 
rate;  whereas,  underestimating  Ropt  is  not  nearly  as  critical.   A  brief 
comparison  between  the  theoretical  and  empirical  results  is  as  follows: 
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jJ        Theoretical  Ropt  Empirical  Ropt 

Re=100  Re=50    Re=10  Re=l 

1  1.8574  1.650  1.725    1.625  

2  1.8916  1.875  1.875    1.875  1.700 

3  1.9064  1.850  1.850    1.825  1.800 

From  the  above  summary,  it  can  be  seen  that  the  theoretical  value  for 
Ropt  is  usually  greater  than  the  empirical  value.   This  situation  might 
lead  to  a  poor  choice  for  Ropt;  therefore,  trying  to  apply  the  theoreti- 
cal results  from  simple  problems  to  more  complex  one  might  lead  to  a 
bad  choice  for  Ropt. 

D.   CONSERVATIVE  AND  TRANSPORTATIVE  PROPERTIES 

One  property  which  is  desirable  in  a  finite  difference  method  is 
the  conservative  property.   A  finite  difference  scheme  possesses  the 
conservative  property  if  it  preserves  a  certain  integral  relation  of 
the  momentum  equations  [Ref.  7].   This  integral  relation  states  that 
over  some  region  the  time  rate  of  change  of  a  given  quantity  must  equal 
the  net  flux  of  the  quantity  across  the  boundary  plus  the  production  rate 
of  the  quantity  within  the  region.   Whether  a  finite  difference  scheme 
preserves  the  conservative  property  depends  upon  the  form  of  the  momentum 
equations  used  and  the  numerical  scheme  employed.   The  possession  of 
this  property  does  not  imply  that  the  scheme  is  more  accurate  than  one 
that  does  not  possesses  it  unless  one  of  the  criteria  of  the  problem 
is  the  conservation  of  a  given  property.   Upwind  differencing,  used  in 
the  simple  cavity  flow  problem,  conserves  vorticity.   The  Arakawa 
Jacobian  conserves  vorticity,  vorticity  squared,  linear  momentum  and 
kinetic  energy;  however,  all  procedures  that  employ  the  Arakawa  Jacobian 
do  not  necessarily  conserve  these  properties. 
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Figure  2.2   Plot  of  overrelaxation  factor  vs.  number  of  passes  required 
for  convergence  for  Re  =  1,  10,  50,  and  100  and  aspect  ratio 
of  1. 
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Figure  2.3  Plot  of  overrelaxation  factor  vs.  number  of  passes  required 
for  convergence  for  Re  =  1,  10,  50,  and  100  and  aspect  ratio 
of  2. 
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Figure  2.4  Plot  of  overrelaxation  factor  vs.  number  of  passes  required 
for  convergence  for  Re  =  1,  10,  50,  and  100  and  aspect  ratio 
of  3. 
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Another  desirable  property  is  the  transportative  property.   The 
finite  difference  form  of  a  flow  equation  possesses  the  transportative 
property  if  the  effect  of  a  perturbation  on  a  transport  property  is 
advected  only  in  the  direction  of  the  flow  [Ref.  7],   An  intuitive 
approach  to  the  problem  leads  one  to  agree  that  a  perturbation  in  the 
flow  field  should  only  move  in  the  direction  of  the  flow.   The  method 
of  upwind  differencing  possesses  the  transportative  property. 

E.   CONVERGENCE  AND  STABILITY 

The  mathematical  basis  for  stability  and  convergence  lies  in  .the 
theory  of  linear  partial  differential  equations.   Linear  theory  is 
used  as  a  guideline  for  analyzing  non-linear  equations  such  as  those 
used  in  this  study. 

A  convergent  finite  difference  scheme  is  one  in  which  all  values  of 
the  finite  difference  solution  approach  the  values  of  the  exact  solution 
as  the  finite  difference  size  approaches  zero.   All  the  terms  in  a  finite 
difference  representation  can  approach  the  corresponding  terms  in  the 
analytic  equation,  and  yet,  it  is  possible  that  the  entire  equation  will 
not  approach  the  exact  solution;  therefore,  the  convergence  criterion 
will  be  unfulfilled.   Convergence,  therefore,  is  concerned  with  the 
limit  of  the  entire  equation  and  not  the  individual  terms. 

Stability  is  achieved  when  the  cumulative  effect  of  all  round-off 
errors  is  negligible.   This  implies  that  the  errors  do  not  increase 
exponentially.   In  the  time  dependent  approach  used  in  both  the  simple 
cavity  flow  problem  and  the  refined  model,  the  stability  of  the  numerical 
method  employed  to  solve  the  vorticity  equation  depends  upon  the  size  of 
the  time  increment  and  the  degree  of  convergence  of  the  stream  function 
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equation.   A  critical  time  step  can  be  determined  for  each  problem. 
Exceeding  this  time  step  causes  the  equations  to  become  unstable  and 
diverge.   It  is  found  that  the  critical  time  increment  depends  upon  the 
grid  spacing,  Ax  and  Ay,  the  Reynolds'  number  and  the  depth  (in  the 
refined  model) .   The  most  stable  solutions  for  the  simple  cavity  flow 
case  occur  at  Reynolds'  numbers  of  approximately  one  hundred  and  stability 
decreases  as  the  grid  size  ratio  of  Ax  to  Ay  increases.   For  the  refined 
model,  the  greatest  stability  is  found  around  Reynolds'  numbers  of  one- 
tenth.   Two  is  the  only  aspect  ratio  that  is  used  in  the  refined  model. 
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Ill .   GOVERNING  EQUATIONS 

A.   BASIC  EQUATIONS  OF  MOTION 

The  governing  equations  in  this  study  are  derived  from  the  hori- 
zontal momentum  equations  and  the  continuity  equation.   Developments 
of  the  momentum  and  continuity  equations  can  be  found  in  many  books  on 
Physical  Oceanography  such  as:   The  Principles  of  Physical  Oceanography 
[Ref.  5].   Applying  Newton's  second  law  to  a  continuous  volume  of  fluid, 
the  horizontal  momentum  equations  can  be  derived.   The  equation  of  con- 
tinuity can  be  derived  from  the  principle  of  the  conservation  of  mass. 
Starting  with  these  equations,  the  complete  development  of  the  governing 
equations  used  in  this  thesis  is  derived  below. 

1 .   Horizontal  Momentum  Equations 

The  Eulerian  representation  of  the  momentum  equations  desired 
in  this  study  is  derived  from  Newton's  second  law  utilizing  the  following 
simplifying  assumptions: 

a.  The  curvature  of  the  earth  is  negligible  for  the  distances 
considered  in  this  study;  therefore,  cartesian  coordinates  can  be  used 
for  the  coordinate  system. 

b.  Fluid  is  homogeneous. 

c.  Fluid  is  incompressible. 

d.  Hydrostatic  approximation. 

e.  Vertical  component  of  Coriolis  force  is  negligible. 

f.  Turbulent  stresses  are  proportional  to  the  gradient  of 
the  mean  flow. 

These  assumptions  imply  that  the  only  external  forces  acting  are:   the 
Coriolis  force,  gravity,  wind  stress,  and  bottom  stress.   Assuming  that 
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primes  denote  non-integrated  variables,  the  Horizontal  Momentum 
Equations  are 

9u'  m      t  3"'  .   .  3uf    f    ,     1  3P'    A  „2  ,    ,  92u'  ,.  1N 

vr—  +  u'  t —  +  v' fv'  = —  +  A  7  u'  +  A  — -  (3.1) 

3t       3x       3y  p  3x     H        v,,  2 

o  Z 


3v*     ,  3v'     ,  3v'    .  ,      1  3P1    .  _2  ,    ,   3"v' 

T—  +  u1  - —  +  v'  - —  +  fu'  = | —  +  A  V  v'  +  A  — —  (3.2) 

3t       3x       3y  p  3y     H        v  „^2 


3z: 


where  in  Equation  3.1 
3u' 


,  3u'     ,  3u' 
u'  - —  +  v'  - — 
3y 


u              - 

3x 

fv' 

1    3P' 

p    3x 

vv 

A       9    U 

v      ^    2 

=  Local  rate  of  change  of  velocity  with  time 

=  Non-linear  field  accelerations 

=  Coriolis  term 

=  Pressure  gradient  term 

=  Lateral  shear  stress 
=  Vertical  shear  stress 


3z 


2 .   Continuity  Equation 

The  Equation  of  Continuity  is  needed  as  one  equation  in  a 
system  of  three  equations  in  three  unknowns.   Physically,  this  equation 
imposes  the  condition  that  no  two  fluid  particles  can  occupy  the  same 
space  at  the  same  time.   Continuing  the  convention  that  primes  denote 
non- integrated  variables  and  neglecting  any  changes  in  the  vertical,  the 
Equation  of  Continuity  can  be  expressed  mathematically  as 

|e.  +  leu:  +  saL  .  0  0.3) 

3t    3x     9y 

For  incompressible,  homogeneous  flow,  the  density,  p,  is  constant; 
therefore,  Equation  3.3  becomes: 
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|^  +  |^=0  (3.A) 

dx    9y 

3.   Integrated  Equation  of  Motion 

Equations  3.1,  3.2,  and  3.4  represent  three  equations  in  three 
unknowns:   u,  v,  and  p.   The  solution  to  these  equations  is  fully  deter- 
mined when  proper  boundary  and  initial  conditions  are  specified.   Two 

approaches  may  be  used  in  order  to  solve  these  equations  for  the  desired 

1 
xntegrated  form. 

a.  The  equations  can  be  integrated  over  the  entire  water 
column,  but  this  cannot  be  done  without  the  non-linear  terms  being 
known  explicitly  as  a  function  of  depth. 

b.  The  real  problem  can  be  replaced  by  one  which  assumes  a 
homogeneous  fluid  and  the  last  term  in  Equations  3.1  and  3.2  can  be 
replaced  by  a  body  force  in  which  u  and  v  are  independent  of  depth. 

The  second  method  is  employed;  therefore,  the  last  term  in  both 
equations  is  replaced  by  a  body  force  which  results  from  the  wind  stress 
and  bottom  stress,  defined  as 

wx    bx      .2  » 

m zr~  =  A  — o  (3.5) 

ph    ph     v  8z2 

wy    by      „2  ? 

— —  -  A (3.6) 

ph    ph     v  3z2 

Substituting  these  terms  into  the  horizontal  momentum  equations  and 
dropping  the  primes  to  denote  that  the  variables  are  now  integrated  over 
the  entire  water  column,  the  integrated  equations  of  motion  are 

9u  ,    9u     9u    _      1  9P  ,  .  _2    tWX  t  ,-  7v 

— -  +  u  -r—  +  v fv= T~  +  A_TV  u  +  — r—         (3.7) 

3t     9x     9y  p  9x    H      pn   ph 


The  following  discussion  is  taken  from  Ref.  4,  pp  12. 
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3v  ,    3v  ,    3v    .      19P    A   2    t"7   Tby 

-T-+U  —  +  V  —  +  fu  = —  +  A  7  v  +  — —         (3.8) 

3t     3x     3y  p  3y    H      ph    ph 

Integrating  the  equation  of  continuity  over  depth,  and  applying 

Leibnitz  Rule  of  integration  yields  the  desired  integrated  continuity 

equation.   If  the  sea  surface  is  assumed  to  be  of  constant  height, 

arbitrarily  selected  as  0,  and  if  h  equals  the  depth  of  the  ocean,  the 

equations  take  the  following  form 
o  o 


J£4"  +  /#*»-° 


-h         -h 
Using  Leibnitz  Rule  where  u  and  v  are  independent  of  depth 


o  o 


|-  f  udz  +  |- 


/  udZ  +  ly"  /  UC 
-h         -h 


The  desired  equation  is 


3  (uh)  +  jLMO  .  0 

3x       3y 


B.   DEVELOPMENT  OF  VORTICITY  EQUATION 

1.   Derivation 

Starting  with  the  equations  of  motion  (3.7  and  3.8)  and  the 

integrated  continuity  equation  (3.9),  the  desired  vorticity  equation 

with  the  wind  stress  assumed  negligible  is  derived  below.   The  local 

wind  stress  is  assumed  negligible  as  this  force  is  not  being  considered 

because  the  oceanic  current  is  the  only  forcing  function  being  examined 

in  this  study.   Assuming  the  wind  stress  is  negligible  implies 

wx    wy 
t   =  t   =  0 

Cross  differentiating  the  equations  of  motion, 

3(3.8)  _  3(3.7) 
3x       3y 
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yields : 


3       ,3v        3u.     ,    3u    ,3v        3u.  3       ,3v        3u,         3v    ,3V        3u. 

3t      3x        3y  3x      3x        3y  3x      3x        3y  3y      3x        3y 

,         3       ,3v        3us     ,     r   3u  3f  3v    ,         3f  /0   ,_. 

+   V  —    (t — )    +    f  —  +   u  ^—  +    f  -r-  +   V  T—  -  (3.10) 

dy      dx        3y  3x  3x  3y  3y 

a    v-    t^y.        ^u\         r    3    rJL—A        _ <L   el ^  i 

V      C3x        3y;         l3x    Qph    '        3y    Cph   )] 

The  vortlcity,    £,    is   defined   as 

C  =  —  -  —  (3.11) 

^        3x        3y  ^'    ^ 

Substituting  into  Equation  3.10  and  collecting  terms  gives 


(3.12) 


Expanding  the  continuity  equation  produces 

3uh   3vh    /3u  ,3V,,,    3h  ,    3h 

-z +  -x =  ("T-  +  -r-)  h  +  UT—  H-  v  —  =  0 

3x    3y     3x   3y        3x     3y 

which  in  turn  yields: 

3u        3v  1    ,      3h  3tu  ,      . _. 

-5—  +  t—  =-r-(uT—  +v  — )  (3.13) 

3x        3y  h  3x  3y 

Substituting   into  Equation   3.12   it   is   found 

If  + » h  »  +  f> +  v  I7  <?  +  f)  -  £  <u  i +  v  t>  (?  +  f)  ■ 

V  ?  +  F  <fc*  "  fe  ^  (3-14) 

now. 
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a    ,t+f\  3      ,S+fN        1        3      /  .n    ,    1        3      /      ,    ^ 

u  ai  (ir}  +  v  17  (— > =  h  u  ai <c  H' f)  +  ir v  17  (c  +  f) 

-^•(u^  +  v"87)  a  +  f)  (3-15) 

therefore: 

9?     L  u    r       9       ^+^  S       ,C+fN1  .    „2       .     3       ,xbXN       3       ,xby 

Tf-  +  h    [u  —   (tt—)    +  v  —    (■*—-)]    =   A.V    ?   +  —    (-r-)-  T~    ("IT") 
3t  3x       h  3y        h  II  3y      ph  3x     ph 

(3.16) 

Introducing  the  stream  function,  ^  ,  and  defining  its  partial  derivatives 

as : 

3y 

3*    h  (3-17) 

^=  Vh 

they  can  be  substituted  into  the  vorticity  equation  to  give 

f  "  F  f  O  +  f  f  ^  -  V*'  +  I"  (V  -  f  c£>   (3-18) 

3t    3y  3x   h     3x  3y   h      H      3y   ph     3x  ph 

If  the  Jacobian,  J,  is  denoted  as 

„,.    „N    3A  3B    3A  3B  ,„  ,„. 

J(A,B)  =  r-  —  -  —  —  (3.19) 

3x  3y   3y  3x 

Equation  3.18  can  be  written 

H  +  j  <♦.  <**»  -  V2^ + 17  (f  >  -  k  (f  >  <3-20) 

2.   Nondimensionalization 

The  equations  are  nondimensionalized  to  generalize  the  results, 
freeing  them  from  the  constraints  of  dimensionality.   Nondimensional 
equations  can  be  written  in  many  different  forms  depending  on  the  way 
in  which  the  nondimensional  parameters  are  formed.   For  example,  non- 
dimensional  time  can  be  defined  as: 
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I-  l 


L/U 


or 


T  = 


l2/Ah 


where 

L  =  Characteristic  length 

U  =  Characteristic  velocity 

A^  =  Horizontal  eddy  diffusivity 
The  choice  of  methods  depends  upon  the  desired  result  which  in  turn 
depends  upon  the  way  in  which  the  equation  is  going  to  be  solved. 
Physical  interpretations  of  the  parameters  developed  may  also  dictate 
the  method  to  be  used.   This  is  because  certain  parameters  may  be 
desired  which  can  be  used  to  help  interpret  the  results. 

In  nondimensionalizing,  reference  quantities  must  be  defined. 
The  following  reference  quantities  are  the  ones  used  in  this  study: 

L  =  Characteristic  length 

U  =  Characteristic  velocity 

oo  J 

D   =  Characteristic  depth 
also  used  in  the  development  are: 

0  =  Beta-plane  approximation  to  the  Coriolis  force 

A^  =  Horizontal  eddy  diffusivity 
Using  these  quantities  the  following  relationships  between  nondimensional 
and  dimensional  quantities  are  derived 

Length  scale 

x-r         '-? 

Depth  scale 

H  -  ± 
H  "  D 
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Velocity  scale 

U  -2-      V  =2- 

u        u 

00  o 

Time  scale 


T-  ' 


L/U 


Nondimensional  Coriolis  parameter 

The  remaining  quantities  to  be  nondimensionalized  can  be  done  so  in  terms 
of  the  above  nondimensional  quantities.   They  are 

Nondimensional  partials  of  the  stream  function 
Knowing  that 

|t  =  -  uh  =  -  UU  HD 

3y  °° 

then 

II  -  r   X  \   M 

9Y  =  lU  D;  3y 


Nondimensional  vorticity 

From  the  definition  of  the  vorticity 


_  _3v  _  _3u  _  .  _a_V   3U,  j» 
C  "  8x   9y  =  l3X   3YJ  L 


thus 


z  -  $->  c 


Nondimensional  bottom  stress 

From  Appendix  A  the  following  relationship  is  obtained 


pL   h  9y   HD   ^U°°U;  3Y 
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therefore 


b     -     b 

—  =  ( J— )  — 

pH    V3LU  '  ph 


Nondimensional  Del  operator 

The  operator  is  defined  as 
2     2 

2     2 

3x    3y 


thus 


2     2   2 

VZ  =  (L  )  VZ 

Notation  is  changed  at  this  point  for  clarity  of  presentation.   Sub- 
scripts are  now  used  to  denote  differentiation.   This  notation  transforms 
Equation  3.20  into 

r+f  r+f         ?  Tbx       xbx 

cT  -*  <V>  +  *x  <*&  ■  v  ?  +  <^r>  - (^ 

x  y  y        x 


All  of  the  nondimensional  terms  are  now  substituted  into  the  dimensional 
1 


equation. 
U 


zt  (r)2 "  U~D^ 


(U  /L)Z+3LF 


DH 


-  +  U  DH'x 
x  L    °° 


(U  /L)Z+BLF 


DH 


1 

y  l 


«>  2 

*h  3  v  z  +  T 


6LU     bx       by 


factoring  yields 

U   2   U   2 

00  CO 


Z+%-F 


H 


x 


U   2 


Z+^F 


H 


bx 


by 


AHU°°   2 

-V v  z  +  eu  i  (V-)  -  <V->  i 

T  3  °°   pH  y    pH  x 


Y  = 
x 


dividing  by: 


U  2 

CO 


For  notational  purposes,  small  x's  and  y's  will  be  used  to  denote  non- 
dimensional  space  coordinates. 
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and  defining  the  two  nondiraensional  parameters 

LU 

00 

Re  =  — —  =  Reynolds  number 
e  =    ■  =  Reciprocal  of  Rossby  number 

CO 

the  equation  becomes 

74-cF  7+pF  1         ?  xbx  rby 

zt  -  *y  «*¥>, +  \  'f1,  ■  h v  z  +  E  [$r>,  -  §r>x]     (3-21) 

Employing  the  definition  of  the  Jacobian, 

7+rF     1    ?        rbx      rby 
ZT  +  J  (I,  (*£££»  -  ^  V2Z  +  e  [(^r)y  -  (^-)x]  (3.22) 

This  form  is  compact  and  neat  for  notational  purposes  but  it  is  not 
suitable  for  solving  numerically  because,  it  is  hard  to  separately 
control  the  individual  forces.   Expanding  the  Jacobian  transforms  the 
equation  into  one  which  is  more  suitable.   Expansion  of  the  Jacobian 
yields 

j  /¥  (Z+e¥))    =  w   fZ+eF)   -  w      (Z+cY) 

J  ^'  lH   jj     x  K   H   ;y    y  K   H   ;X 

H(Z+eF)  -(Z+eF)H         H(Z+eF)  -(Z+eF)H 
=  y   r 1 Y_i  _  y   r * X-i 

x  l        h2         ]     y  [        h2         ] 

Y  (Z+eF)     Y  (Z+eF)     ,7J_  _. 

H  H         7        L  y  x    x  yJ 

rl 

(f  2  -f  Z  ) 

=        X  y„  y  X     +  £  (Y   F  -¥   F   )   +    (Z+^F)    [H  !  -H  ?   ] 
H  Hxyyx  u2  xyyx 

rl 

4J  C^Z)  +  f  J  (V,F)  +  C^^)  J  (H,f) 
H  H  R2 

Substituting  this  form  into  Equation  3.22 
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ZT  +  I  J  (¥,Z)  +  f  J  (f  ,F)  +  (^)  J  (H,V)  = 

H. 

(3.23) 

,    9         bx       by 

vzz  +  e  [(ttt-)„  -  Cr5-)J 


Re  pH  y    pH  x" 


where 


Z  =  Local  rate  of  change  of  vorticity  with  time 

—  J  (Y,Z)  =  Advective  term 
H 

—  J  (y,F)  =   Planetary  vorticity  tendency 
H 

(Z+eF) 


J  (H,H0      =  Topographic  vorticity  tendency 
H 

1   2 

—  V  Z  =   Frictional  term 

Re 

bx       by 

T  T 

e  [  (~rr~)      -    (~^~)     I  =  Bottom  friction 
pH  y    pH  x 

This  form  is  much  more  convenient  for  deleting,  weighting  or  evaluating 

the  acting  forces  individually.   The  method  for  incorporating  the  bottom 

friction  will  be  explained  in  the  section  on  the  refined  model  and  the 

complete  development  of  the  bottom  friction  is  contained  in  Appendix  A. 

3.   Nondimensional  Parameters 

Two  nondimensional  parameters  occur   in  the  vorticity  equation. 

The  first  one,  the  Reynolds'  number,  is  defined  as 

LU 

CO 

where  L,  U  ,  and  A^  have  been  defined  previously.   Physically  the 
Reynolds'  number  is  a  ratio  of  the  viscous  to  the  inertial  forces. 
The  second  parameter,  e,  is  given  by 

*"    U 
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This  can  be  interpreted  as  the  reciprocal  of  the  Rossby  number  which  is 
a  measure  of  the  effect  of  the  Coriolis  force. 

The  solution  of  the  equations  depends  only  on  the  values  of 
the  parameters,  not  on  the  values  of  the  individual  terms  in  the  param- 
eters.  Nevertheless,  estimates  must  be  made  of  the  individual  terms  in 
order  to  determine  the  range  of  values  for  the  nondimensional  parameter 
which  should  be  investigated.   Certain  terms,  in  particular  3  and  to  some 

degree  A^,  are  fairly  well  known.   The  horizontal  eddy  diffusivity,  iL  , 

7       9    2 
has  been  found  for  turbulent  flow  to  be  between  10   and  10   cm  /sec. 

7   2 
The  value  of  10   cm  /sec  will  be  used  in  this  development.   $  can  be 

developed  from  its  definition. 

2fi  cos  9 
o 


where 

fi  =  Angular  rotation  of  the  earth  =  7.29  x  10   rad/sec 
9  =  Latitude  =  37   (for  this  study) 

o 

r  =  Radius  of  earth  =  6 . 38  x  10   cm 
hence, 

a   -   2(7.29  x  10"5)  x  .8   .  p_   _  -13   -1    -1 
P  =  ~ =  1.83  x  10    cm   sec 

6.38  x  10 
Determining  representative  values  for  L  and  U  which  are  the  nondimensional 

oo 

length  and  velocity  scales  is  not  as  easy.   Physical  measurements  can  be 
relied  upon  to  give  an  estimate  for  U  which  can  be  considered  the  free 

oo 

stream  velocity  or  in  other  words  the  average  velocity  of  the  current 
along  the  coast  which  is  unaffected  by  the  presence  of  the  cavity.   The 
best  estimate  that  can  be  made  indicates  that  a  value  of  approximately 
one-tenth  to  one  centimeter/second  is  a  reasonable  one.   Theoretical 
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considerations  are  relied  upon  to  obtain  a  representative  value  for  L. 
The  value  chosen  for  L  is  the  most  questionable  of  all  the  character- 
istic values.   A  value  of  one  mile  (1.61  x  10   cm)  is  chosen  since  it 
represents  the  smallest  possible  size  for  any  feature  to  be  distinguish- 
able on  the  grid  employed.   This  is  because  one  mile  is  the  size  of  one 
grid  space  on  the  arrays  used  to  solve  the  equations  numerically.   Using 
these  values,  it  is  found  that  approximate  values  for  the  nondimensional 
parameters  used  in  this  study  are 

Re  =  1.61  x  10~2 

e  =  4.74  x  10"3 

C.   DEVELOPMENT  OF  STREAM  FUNCTION  EQUATION 

The  vorticity  equation  is  one  equation  in  two  unknowns;  the  vorticity, 
C,  and  the  stream  function,  ty.      In  order  to  solve  the  vorticity  equation, 
another  equation  relating  to  one  or  both  of  these  variables  must  be 
found.   A  second  equation  is  derived  from  the  definitions  of  the 
vorticity  and  stream  function  and  equating  the  two. 

1.   Derivation 

Respectively,  the  vorticity  and  the  derivatives  of  the  stream 

function  were  previously  defined  as 

C   =  v     -   u 
x  y 

il>     =  vh 
x 

\b      =   -uh 

y 


therefore 


x  y 


<  -  <-f  >  +  H> 
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thus  , 


h\b        -   ib  h  hil;        -   i>  h 

_        xx x  x 


2  2 

h  h 

finally,   separating   terras 

C  =  £  V2^  -  —  (if,  h     +  1(1  h  )  (3.2*) 

h  ,  2        x  x  y  y 

Equation  3.24  can  be  solved  for  the  stream  function  if  it  is  converted 
to  its  numerical  representation  by  substituting  the  proper  relationships 
for  the  derivatives.   This  process  is  carried  out  in  the  numerical 
approach  sections  which  occur  in  both  the  simple  cavity  flow  problem 
and  the  refined  model. 

2.   Nondimensionalization 

Following  the  same  procedure  as  in  Section  III-B-3  and  using  the 
same  nondimensional  terms  as  used  for  the  vorticity  equation,  the  stream 
function  equation  can  be  nondimensionalized.   Substituting  into  the 
stream  function  equation  yields 

Z  (J=)  -  V2  £->  ,  (U.LD)  |  (i)  -  Jj  (ij)  [Tx(UJ»  Hx  <2)  + 

L  n    D 

\    <U~D)    Hy    <ff>] 
collecting   terms 

U  „2UJ      ULD  .  U   D2 

CO  V      Y  oo  oo 

Z     (— )     =    ^    (-=— )     "   ^    [?    H       +    ?    H     ]     (— — ) 

L      H   L2D     H2   x  x    y  y   L])2 
therefore,  all  characteristic  values  cancel  giving 

Z  =  ^  V2¥  -  =rr  (V  H  +  Y  H  )  (3.25) 

H       R2   x  x    y  y 

This  is  the  desired  nondimensional  stream  function  equation  which  is  used 
in  conjunction  with  the  vorticity  equation  to  solve  for  both  the  stream 
function  and  the  vorticity. 
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IV.   SIMPLE  CAVITY  FLOW 

A.   GOVERNING  EQUATIONS 

The  governing  equations  for  the  simple  cavity  flow  problem  derive 
directly  from  the  equations  developed  in  Section  III.   Two  equations 
are  necessary:   The  vorticity  transport  equation,  Equation  3.23,  and 
the  stream  function  equation,  Equation  3.25. 

ZT  +  |  J  (¥,Z)  +  |  J  <Y,F)  +  IS+SEl  J  (H,Y)  = 

H 

Z  =  1  v2^  -  -^  (V  H  +  ¥  H  ) 
H       R2   x  x    y  y 

The  following  assumptions  are  made  in  the  simple  cavity  flow  problem 

1.  Bottom  friction  is  negligible 

2.  Depth  is  constant 

3.  Coriolis  force  is  negligible 

Using  these  assumptions,  Equations  3.23  and  3.25  become 

ZT  =  Re"  V^  "  J  (VF'Z)  (4,1) 

Z  =  V2Y  (A. 2) 

where  the  stream  function  for  a  constant  depth  reduces  to 
V  =  -U 

y 
¥  =  v 

x 

In  developing  the  general  equation,  it  is  found  that  for  convenience 

of  form  it  is  better  to  define  the  vorticity  as 

Z  =  V  -  U 
x    y 
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although  the  normal  convention  for  the  vorticity  is 

Z  =  U  -  V 
y    x 

Simplifying,  the  equation  eliminates  the  need  for  this  unconventional 

form;  therefore,  the  normal  convention  for  the  vorticity  is  used. 

Changing  the  convention  changes  the  definition  for  the  velocity  in  terms 

of  the  stream  function  and  some  of  the  signs  in  the  vorticity  equation. 

The  partials  are  now  defined  as 

y  =  U 

y 

f  =  -v 

X 

Using  this  convention  and  the  fact  that  J(A,B)  =  -J(B,A)  the  vorticity 
equation  becomes 

ZT  =  Re"  V^  "  J  (ZsH,)  (4,3) 

The  above  equation  for  the  vorticity  is  in  suitable  form  for  solving 

by  utilizing  the  Arakawa  Jacobian,  but  a  slightly  different  form  is 

needed  to  solve  by  upwind  differencing.   Applying  the  definition  of  the 

Jacobian  we  have 

J(Z,y)  =Z!   -  Z  ¥ 
x  y    y  x 

substituting  for  the  partials  of  the  stream  function 

J(Z.Y)  -  Z  U  +  Z  V 

x     y 

this  term  can  be  expanded  to 

UZ  +  VZ   =  (UZ)   +  (VZ)   -  Z[U  +  V  ] 
x     y       x       y      x    y 

but  the  bracketed  term  equals  zero  by  the  continuity  equation,  and  it  is 

found  that  the  Jacobian  becomes 

J(Z,Y)  =  (UZ)   +  (VZ)  (A. 4) 

x       y 
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Introducing  the  new  term  into  the  vorticity  equation,  the  desired  form 
is  obtained. 

ZT  =  Re"  y2z  "  [(UZ)x  +  (VZ)y]  (4,5) 

B.   NUMERICAL  APPROACH 

1 .   Upwind  Differencing  for  Solving  the   Vorticity  Equation 

Upwind  differencing  achieves  static  stability  by  differencing 
in  the  direction  of  the  flow.   Differencing  of  the  advective  terms  in 
this  method  is  always  upwind  of  the  point  in  question.   Considering 
Equation  4.5,  forward  time  differencing  can  be  used  for  the  first  term 
and  the  technique  for  representing  the  Laplacian  can  be  used  for  the 
second.   In  is  the  advective  terms  in  brackets  which  must  be  handled 
in  a  special  manner. 

Zn+!  -  Zn     ,    Zk    .+Zk+|  .-2Zk     Zk  _+Zk+1   -2Zk 

(Ax)  (Ay) 

-  [advective  terms]  (A. 6) 

where 

n  is  time  step  level 
k  is  iteration  level  in  space 
The  numerical  scheme  used  to  represent  the  advective  terms  depends  upon 
the  direction  of  the  flow,  sign  of  U  and  V,  and  on  whether  or  not  the 
flow  reverses  itself  between  two  successive  grid  points.   Using  upwind 
differencing,  the  direction  of  the  flow  and  occurrence  of  flow 
reversals  must  be  checked  at  each  grid  point.   One  upwind  differencing 
scheme  is  as  follows : 
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NO  FLOW  REVERSALS  OCCUR 


+  U 


+  V 


-  U 


-  V 


(uz).  -  (uz)i_1 

Ax" 

(vz).  -  (vz) 

Ay 
(uz),+1  -  (UZ). 

Ax 
(uz)j+1  -  (UZ)j 

Ay 


Backward  differencing 
Backward  differencing 
Forward  differencing 
Forward  differencing 


FLOW  REVERSALS 

There  are  several  methods  which  adequately  handle  flow  reversals 
The  one  incorporated  into  this  study  is  an  averaging  scheme.   Consider- 
ing the  x-direction  advective  term,  it  is  represented  by 

(URZR  "  W 


Ax 


where 


UR  = 


Ui+1  +  Ul 


Z  =  Z   if  U  >  0 
R    i     R 


\  '   Zi+1  "  UR  <  ° 


U.  +  u.  . 


Z,  -  Z.   'if  U  >  0 

L    l-l     L 


ZT  =  Z.  if  UT  <  0 
L    l     L 


The  formula  for  the  y-direction  term  is  analogous  to  the  x-direction 

formula. 

Interpreting  the  situation  physically,  a  flow  reversal  appears 

as  an  artificial  source  or  sink  depending  on  the  sign  (-  to  +  is  a 

source  of  vorticity) .   Averaging  eliminates  the  source  or  sink. 

2.   Method  for  Solving  Vorticity  Equation  Employing  the 
Arakawa  Jacob ian 

This  approach  to  the  problem  is  very  simple.   Forward  time 

differencing  is  used  for  local  rate  of  change  of  vorticity  with  time, 
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the  representation  for  the  Laplacian  is  used  for  the  frictional  term  and 
the  representation  for  the  Jacobian  is  the  one  developed  by  Arakawa. 
Originally  only  upwind  differencing  was  going  to  be  used  to  solve  the 
simple  cavity  flow  problem  but  when  the  aspect  ratio  was  increased  to 
three  the  upwind  differencing  method  would  not  converge.   Utilizing 
the  Arakava  method,  solutions  were  obtained  but  at  the  expense  of 
decreased  detail.   The  following  is  the  numerical  equation  which  was 
derived  for  this  procedure. 


zn+1  -  zn 


At 


^  =  ir  [ 


k                  k+1  k 

Zr,,     .   +  ZK  ,     .    -   2Z. 
1+1  >  J 1"1?J Lj 


k+1  k 

+  Z.     .    ,    -   2Z. 


Re 


(Ax) 


1  +     i,j+l  i,.j-l 

(Ay)2 


i*l]       (4.7) 


12AxAy 


1+1,3  i,J-l      i+l,j~l      x,j+l      1+1,  j+1 

V.     .,-  (Z,,.     ,+Z.,,     .,,-Z.    .     .-Z.    _     ,,_)   + 

i,j+l  i+l,3      1+1,3+1      i-l,3      1-1,3+1 

Y.       .  (Z.    .,,+z.        ..,-z,    .   ,-z,   .   ,  .)  + 

1-1,3  1,3+1      1-1,3+1      1,3-1      1-1,3-1 

v.    .   ,  (z.    ,    ,+Z.        .   --Z..-    ,-z,.-    .      )  + 

i,j-l  i-l,3      1-1,3-1      i+l,3      1+1,3-1 

Vl,j+1    (Zi41,j-Zi,j+l)   + 


v.  ,    ,.,    (z.    ..,-z,   .    .)  + 

1-1,3+1        i,j+l      i-l,3 


y4  ,  .  ,  (z.  1  ,-z,  .  j  + 

1-1,3-1        1-1,3      1,3-1 


Yi+i,j-i  (zi,j-rzi+i,j} 


3.   Optimum  Overrelaxation  to  Solve  Poisson's  Equation  for  the 
Stream  Function 


The  numerical  equation  used  to  solve  for  the  stream  function  is 
derived  in  the  following  manner. 

Z  =  V2f 
Using  the  representation  for  Laplacian 


n      =     i+1>J         j-kii  i»i  +    i»J+1        i>J"1  i'i 


i,j  (Ax)2  (Ay) 


Letting  B  =  Ax/Ay 

(Ax)2  Zn  ,  =  ^kin  ,  +  «^+J  .  -  2Tk+l  +  B2  [Yk  .,,  +Vk+]   ,  -  24>k+1] 
this  yields 

vk+l ~r-  rf,i    .+^   .+B2(vFk  .11+vfk+1  J-(Ax)2zn   .] 

1,3    2(1+B  )    i+1»J   1_1'J      1,3+1   1,3-1        1,3 
This  method  for  deriving  an  equation  for  the  stream  function  is  known 
as  Richardson's  method.   If  new  values  are  used  as  soon  as  they  are 
available,  it  becomes  Leibmann's  method.   A  further  refinement  is 
achieved  by  using  optimum  overrelaxation.   The  formula  for  this 
technique  is 

yk+1  =  yk   +  Ropt    [H,k    +H,k+1  +B2(lJ/k    +4,k+l   ) 
i,3     i,3    2(1+B2)    i+1'J   i-1'3     i,3+l   i,3"l 


(4.8) 


(Ax)2  Zn   -  2(1+B2)yk  .] 
1,3  i,3 


Ropt  is  the  optimum  overrelaxation  factor  which  in  this  study  has  been 
found  to  be  approximately  1.72  for  the  best  overall  results.   The 
bracketed  term  is  the  residual  and  as  the  solution  converges  to  the 
exact  solution  it  approaches  zero. 
4.   Boundary  Conditions 

There  are  many  solutions  to  any  given  partial  differential 
equation.   Boundary  and  initial  conditions  must  be  specified  in  order 
to  determine  a  unique  solution.   Figure  4.1  identifies  the  boundaries 
which  must  be  specified  in  the  simple  cavity  flow  problem. 
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1. 

Inflow 

2. 

Wall 

3. 

Wall 

4. 

Wall 

5. 

Lid 

6. 

Outflow 

Figure  4.1   Specifications  of  Boundaries 


Each  boundary  and  its  associated  equation  are  discussed  separately  in 
the  following  sections, 
a.   Inflow 

(1)  Velocity  Profile 

To  be  able  to  specify  the  inflow,  a  velocity  boundary 
layer  profile  must  be  defined.   Three  velocity  profiles  are  considered: 
a  polynomial,  logarithmic,  and  hyperbolic  tangent.   It  is  found  that 
there  is  little  difference  in  the  results  obtained  by  using  these 
different  profiles.   The  profiles  are  plotted  for  comparison  in  Figure 
4.2.   The  equations  are  normalized  so  that  the  values  range  between 
0  and  the  maximum  velocity,  1.   Over  this  range  the  polynomial  and  the 
logarithmic  profiles  agree  well  although  the  logarithmic  profile  does 
not  approach  the  maximum  value,  1,  asymptotically  as  desired. 

The  first  profile  considered  is  a  polynomial  derived 
experimentally  for  small  scale  flow  in  Ref.  6.   It  is  expressed 
mathematically  by  the  following  fifth  degree  polynomial: 


U  =  .000493  +  2.408031y  -  2.841989y2  +  2.796410y3 
-  1.904102y4  +  .533821y5 


(4.9) 
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The  second  profile  derives  from  Prandtl's  equation 
[Ref.  8].   According  to  Prandtl's  theory  it  can  be  shown  that  the 
turbulent  shearing  stress  becomes 

where 

K  =  Von  Karman's  constant 

u  =  average  velocity 
Directly  from  Equation  4.10 

du  _  K  /  o  1_ 
dy     V  p  y 

If  it  is  assumed  that  the  shearing  stress,  t  ,  is  constant  and  the 

density,  p,  is  constant  then 


K 
'p 


./ —  =  constant  =  c. 
V  p  1 


therefore 


du  _  fl 

dy  "  y 


and 


u  =  c.  In  y  +  c 


To  normalize  this  profile  it  is  required  that 

u  =  0  at  y  =  0 
u  =  1  at  y  =  1 

thus 

c2  =  0 

Cl  =  1/ln  (e  +  1.0) 

y  =  ey'  +  1.0 

where  e  =  2.7182818 
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The  velocity  equation  becomes 

U  =  In  (e\   1.0)  ln  (£y  +  L-0)  (4al) 

The  last  profile  considered  is  the  hyperbolic  tangent. 
This  profile  is  considered  because  it  is  well  behaved  and  possesses 
two  desirsble  characteristics.   First,  it  approaches  1  asymptotically 
and  secondly  near  the  boundary  it  is,  to  a  close  approximation, 
linear.   In  order  to  make  the  value  of  the  profile  vary  from  0  to  1 
over  the  values  of  y  =  0  to  y  =  1,  the  argument  Try  must  be  used. 
The  resulting  equation  is 

U  =  tan  h  (Try)  (4.12) 

(2)  Stream  Function  Profile 

The  stream  function  profile  is  obtained  by  integrating 
the  velocity  profile  from  the  boundary  to  y.   It  is  found  that  unless- 
a  well-smoothed  stream  function  profile  is  obtained  from  the  velocity 
profile,  instabilities  in  the  vorticity  equation  develop.   For  example, 
the  following  technique  does  not  work  well  in  deriving  the  stream 
function  from  the  velocity  profile.   Given 

central  differencing  yields 

i,.i+i Lozi_  u 

2Ay 
where  in  the  preceding  equation  the  subscript  i  =  1  denotes  the  first 
grid  column  which  is  the  inflow  boundary  (Figure  4.1).   Thus 

\i+i -  'i.j-1 +  2AyU 

Using  this  procedure,  a  profile  like  the  one  in  Figure  4.3  is 
developed. 
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Stream  Function 

Figure  4.3  Unsuitable  boundary  layer  profile 
This  profile  causes  the  vorticity  field  to  alternate  in  sign.   For  this 
reason,  all  three  velocity  profiles  examined  are  integrated  analytically 
in  order  to  obtain  the  stream  function  profiles.   The  three  stream 
function  profiles  are 
POLYNOMIAL 


mn/cn    2.408031  2   2.841989   3 
"P  =  .  00049  3y  H ~ y  -  - y 


2.796410  4   1.904102  5   .533821  6 
+ T y = y  + 7 y 


(4.13) 


LOGARITHMIC 


(ey  +  1.0) 
In  (e  +  1.0) 


[In  (ey  +  1.0)  -  1.0] 


(4.14) 


HYPERBOLIC 


Y  =  —  log  (cos  h  iry) 


(4.15) 


The  three  profiles  are  plotted  in  Figure  4.4.   It  is  apparent  how 
similar  the  three  profiles  are.   The  choice  of  the  hyperbolic  profile 
for  the  one  used  in  this  study  is  based  primarily  on  the  fact  that  it 
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possesses  the  two  theoretically  desirable  properties  previously 
discussed. 

(3)  Vorticity 

The  development  of  the  vorticity  input  equation  is  much 
simpler  than  the  stream  function  equation.   Only  one  assumption  is 
necessary  in  this  development ,  that  is 

f  €  >    -  ° 

3x   dX  i  , 

this  can  be  approximated  by 

_3V,     =  _3_V| 

3x'   .    3x'   . 
1,3      2, j 

Remembering  the  relationship  between  the  velocity  and  stream  function 

then 

3U.      32f, 


8y'l,3 

3y2  i,J 

3  V, 
3x'       . 

= 

2 

■s    2', 

3   4? 

„    2 

l,j      3x  l,j      3x   2,j 
Defining  the  vorticity  as 


3Ui     _  _3Vi 


substitution  yields 

2         2 

31,       91 . 
7     =  +  

^    3yZ  l,j    3xZ  2,j 
Centered  differencing  gives 

Z    =    ,J     ,:]      '-1  I  3,;1   1,:1     }1      (4.13) 
1,j  (Ay)2  (Ax)2 
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b.   Walls 

All  of  the  walls  can  be  handled  in  a  similar  manner.   Two 
assumptions  are  necessary  in  order  to  find  the  value  of  the  stream 
function  at  the  walls. 

First,  it  is  assumed  that  the  no  slip  condition  applies  at 
all  the  walls;  therefore,  at  the  walls 

U  =  V  =  0 

V      =  Y  =  0 

y   x 

hence,       V   =  constant 

Secondly,  assuming  that  the  arbitrary  constant  of  integration  is  zero, 

then  for  all  the  walls 

¥  -  0  (4.14) 

The  vorticity  is  handled  in  a  slightly  different  manner. 
Consider  first  the  walls  parallel  to  the  X  direction,  at  the  wall  it 
is  obvious  that  V  =  0,  therefore 

z    m™L\         =  a2yi 

wall   3y'   1n    .2'   n1 
wall    3y  wall 

Expanding  ¥  into  a  Taylor's  series 

°  J»°  ? 

Y       n^i=X        i  i    +  "77  Ay+^r  ;r  (Ay) 

wall+1     /     wall       isy  2    _    2        .  - 

/  /       wall  9y     wall 

this  gives : 

Z  „  -— SSim  (4.15) 

Wal1     (Ay)2 

U  equals  zero  for  the  walls  parallel  to  the  y  direction,  thus 

=   av,   _  A 

wall     8X1   ...    '21   .. 
wall   3x  wall 

Expanding  again  yields: 
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o 

+  |?|      Ax+^-^4|     (Ax)2 


r- 


wall+1   /wall    ^x1   n.      2  „  2 

wall        9x    in 
wall 


therefore 

wall     /,  \2 
(Ax) 

c.  Downstream  Continuation 

The  stream  function  at  the  outflow  is  determined  by  a  simple 
extrapolation  procedure.   Letting  IL  be  the  outflow  boundary  and  IL-1 
and  IL-2  the  two  preceeding  grid  points,  then  the  desired  relation  is 
given  by 

(J)       +  vj; 

_  YIL-2    *IL 
IL-1        2 

thus 

Y   =  2^       -  H*  (A.  ]  7) 

IL     IL-l,j     IL-2,j  V^.J/; 

The  vorticity  at  the  outflow  is  assumed  to  be  transferred 
out  of  the  grid.   This  yields 

ZIL  "  ZIL-1  <*•«> 

This  can  be  interpreted  as  meaning  that  there  is  no  production  of 
vorticity  between  the  last  two  grid  points. 

d.  Lid 

If  the  Lid  is  considered  to  be  a  frictionless  impermeable 

wall  then  it  is  implied  that  it  is  a  streamline.   A  streamline  is  a  line 

along  which  ¥  is  constant;  since,  the  stream  function  is  assigned  a 

value  at  the  inflow  boundary  then  it  follows  that  all  of  the  other 

values  along  the  Lid  must  have  the  same  value,  or 

Y      =  T  (4.19) 

i,LID    1,LID  v    ' 
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Determining  the  vorticity  at  the  Lid  is  based  upon  two 
assumptions  about  the  velocity  at  the  Lid.   They  are 

f|     =0 

3X  LID 

_aUi   =  3U. 

3y  LID  =  3y  LID-l 

from  these 

Z.  TTT.  -  Z.  TTn  (4.20) 

x,LID     i,LID-l 

This  can  be  interpreted  as  a  linear  extrapolation  of  the  U  velocity 

component  up  to  the  Lid. 

e.   Corners 

It  is  desirable  to  retain  the  effect  of  sharp  corners  in 

order  to  approximate  reality  to  as  closs  a  degree  as  possible.   The 

stream  function  is  no  problem  since  the  value  at  the  corner  is  the  same 

as  the  rest  of  the  wall,  or  zero.   The  approach  to  the  vorticity  is  a 

little  more  complicated  and  it  is  handled  by  using  different  values  for 

the  corner  depending  upon  the  point  being  evaluated.   The  scheme  which 

is  used  is  as  follows 


A    2(Ti,i+i  -  Ti,i) 

2 
(Ay)Z 


•l.J+1 


wall 


•  i+l,j 


Ax 


20'       .  -  v.  .) 

B  =  1+1>J LiJ_ 

2 
(Ax)Z 

If  the  vorticity  at  i,j+l  is  being  evaluated  the  value  A  is  used  for 

Z.  .  when  Z.,,  .  is  being  evaluated  the  value  B  is  used.   This  scheme 
i,J       i+l, 3 

retains  the  effect  of  sharp  corners.   A  similar  method  is  employed  for 

the  other  corner. 
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5 .   Convergence  and  Stability  Criteria 

As  is  discussed  in  Section  II-D,  the  convergence  of  the  vorticity 
transport  equation  depends  upon  the  degree  of  convergence  of  the  stream 
function  equation  and  the  time  step  taken.   The  maximum  allowable  time 
step  in  turn  depends  upon  the  grid  spacing,  the  maximum  velocities  and 
the  Reynolds'  number.   This  critical  time  step  can  be  determined  by 
employing  a  discrete  perturbation,  stability  analysis  [Ref.  7].   Using 
this  method  it  has  been  determined  that  the  critical  time  step  for  the 
upwind  differencing  method  is  given  by 

Atc  = 1 (4.21) 

max  +  max   2_  ,   1       1  ■, 

Ax     Ay    Re  '    .2    ,.  .2J 
(Ax)     (Ay) 

The  time  step  taken  must  be  smaller  than  this  critical  step  if  the 
equations  are  to  converge. 

It  is  found,  by  making  many  test  computer  runs,  that  it  takes 
approximately  100  relaxations  of  the  stream  function  in  order  to  allow 
the  vorticity  equation  to  converge  to  within  acceptable  limits.   The 
number  of  relaxations  is  less  for  the  method  employing  the  Arakawa 
Jacobian  than  for  the  upwind  differencing  method  since  the  former 
method  averages  three  values  for  the  advective  term,  minimizing  the 
growth  of  instabilities. 

The  degree  of  convergence  of  the  equations  is  not  very  strict 
since  this  is  a  preliminary  study,  although  it  is  felt  that  the  conver- 
gence obtained  is  great  enough  so  that  any  further  changes  that  might  occur 
is  less  than  the  accuracy  of  the  plotting  capabilities  of  the  Calcomp 
plotters  used.   In  general,  the  equations  are  relaxed  until  the  residual 
is  one-tenth  of  one  percent  of  its  original  value.   The  original  value 
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is  the  one  determined  by  the  first  guess  field  which  is  obtained  by 
making  the  entire  field  outside  of  the  cavity  equal  to  the  input  boundary 
layer  profile  and  all  the  values  in  the  cavity  equal  to  zero.   Using 
these  criteria  and  a  time  step  equal  to  eight-tenths  of  the  critical 
value,  it  takes  between  twenty  and  forty  minutes  of  computer  time  on 
the  IBM  360/67  to  obtain  the  desired  results  for  each  set  of  input 
parameters,  which  are  the  Reynolds'  number  and  the  aspect  ratio. 
6 .   Computational  Sequence 

Simply  stated,  the  procedure  used  to  solve  the  numerical 
equations  involves  calculating  the  vorticity  transport  equation, 
relaxing  the  stream  function  equation,  taking  a  time  step  and  calcu- 
lating the  new  vorticity  field.   This  procedure  is  repeated  until  the 
equations  are  converged  to  within  acceptable  limits.   A  more  compre- 
hensive list  of  the  steps  involved  in  this  procedure  is  as  follows: 
STEP1  -  Define  all  constants  which  specify  the  problem  and/or 

minimize  computational  time. 
STEP2  -  Specify  the  initial  stream  function  by  either  reading  in  an 
initial  guess  field  or  extrapolating  values  from  the  inflow 
boundary  condition. 
STEP3  -  (Needed  only  for  upwind  differencing  method.)   Specify  the 

initial  velocity  fields  utilizing  the  initial  stream  function 
field. 
STEP4  -  Specify  the  initial  vorticity  field  by  either  reading  in  an 

initial  guess  field  or  determining  the  field  by  using  Poisson's 
equation. 
STEP5  -  Calculate  all  the  boundary  conditions  for  the  vorticity. 
STEP6  -  (Needed  only  for  upwind  differencing  method.)   Calculate  the 
velocity  fields  from  the  stream  function  field. 
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k+1 
STEP7   -  Calculate  Z   "  using  the  vorticity  transport  equation. 

STEP8  -  Check  the  vorticity  field  for  convergence.   Exit  to  STEP11 
if  convergence  criteria  is  met,  continue  if  not. 

STEP9   -  Calculate  the  new  *F  distribution  using  optimum  overrelaxation 
holding  the  vorticity  field  constant. 

STEP10  -  Advance  one  time  step  and  return  to  STEP5. 

STEP11  -  Plot  stream  function  and  vorticity  field  computed  on  Calcomp 
plotters . 
7.   Programming  Optimization 

One  of  the  objects  of  solving  the.  simple  cavity  flow  problem  is 
to  gain  experience  in  programming  problems  of  this  nature.   An  important 
aspect  of  programming  is  being  as  efficient  as  possible  in  order  to 
minimize  the  computer  time  necessary  to  solve  the  problem.   Two  basic 
principles  to  follow  in  order  to  do  this  are:   minimize  operations  and 
structure  the  program  efficiently.   Minimizing  operations  can  be  divided 
into  two  categories.   First,  one  can  minimize  computations  by  manipu- 
lating the  numerical  equations  until  they  are  in  the  most  efficient  form 
to  be  programmed.   This  basically  breaks  down  to  collecting  like  terms 
and  defining  new  constants  (to  be  put  in  the  initialization  part  of  the 
program)  in  order  to  eliminate  unnecessary  calculations.   Secondly,  all 
unnecessary  computations  should  be  eliminated.   An  example  will  demonstrate 
how  this  is  different  from  the  first  category.   The  critical  time  step  in 
the  simple  cavity  flow  problem  involves  the  maximum  components  of  the 
velocity  (Equation  4.21).   This  requires  many  time  consuming  computations 
in  order  to  determine  the  maximum  values  on  the  grid.   It  was  found  that 
after  five  time  increments  the  maximum  values  of  the  velocity  components 
varied  very  little  and,  therefore,  the  time  increment  stabilized.   It  was 
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therefore  unnecessary  to  continue  these  time  consuming  calculations 
to  determine  the  critical  time  step  since  it  remained  virtually  constant. 
The  second  basic  principle  of  structuring  the  program  efficiently 
can  make  a  tremendous  difference  in  the  amount  of  computer  time  used.   A 
summary  of  some  of  the  things  this  author  has  found  helpful  is  contained 
in  the  following  list. 

1.  Avoid  subroutines,  especially  in  DO  loops 

2.  Avoid  function  statements 

3.  Avoid  unnecessary  indexing  and  if  an  index  appears  more 
than  once  in  a  calculation,  equivalence  it  to  a  nonsubscripted  variable 

A.   Avoid  mixed  modes 

5.  Avoid  computed  and  assigned  GO  TO  statements 

6.  Convert  divisions  to  multiplications  where  possible 

7.  Initialize  values  in  DATA  statements 

Most  of  the  things  in  the  list,  which  it  is  advisable  to  avoid,  were 
developed  for  the  convenience  of  the  programmer,  but  the  price  that  is 
paid  for  this  convenience  is  increased  computational  time. 

C.   RESULTS 

The  complete  results  of  the  simple  cavity  flow  problem,  consisting  of 
twenty-four  Calcomp  plots  of  the  vorticity  and  stream  function,  is  con- 
tained in  Appendix  B.   These  plots,  twelve  stream  function  and  twelve 
vorticity,  are  for  Aspect  Ratios  of  1,  2,  3,  and  Reynolds'  numbers  of  1, 
10,  100,  and  1000.   The  plots  are  of  only  a  portion  of  the  total  grid  as 
shown  in  the  following  diagram  in  which  the  shaded  area  is  that  area 
which  is  plotted. 

Contours  are  not  plotted  at  equal  intervals;  therefore,  the  relative 
magnitude  of  the  circulation  within  a  plot  cannot  be  determined  from  the 
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stream  function  gradients.   Only  the  direction  of  the  flow  is  indicated 
by  the  streamlines  since  the  velocity  vectors  are  tangent  to  the  stream- 
lines.  All  plots  are  plotted  using  the  same  contour  levels  so  that  a 
realistic  comparison  between  plots  can  be  made. 
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PORTION  OF  TOTAL  GRID  PLOTTED 


Figure  4.5  Portion  of  numerical  grid  plotted 

The  main  purpose  of  this  part  of  the  investigation  is  an  examination 
of  the  different  circulation  patterns  which  might  occur.  With  this  pur- 
pose in  mind,  four  of  the  more  interesting  flow  patterns  are  presented 
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in  Figure  4.6.   An  indication  of  the  effect  of  increasing  the  magnitude 
of  the  flow  by  increasing  the  Reynolds'  number  is  shown  in  the  top  two 
plots.   Not  only  is  the  magnitude  of  the  flow  increased  when  the  Reynolds' 
number  is  increased,  but  the  pattern  becomes  more  circular  and  the 
center  of  the  gyre  is  positioned  closer  to  the  mouth  of  the  cavity.   The 
lower  diagrams  represent  patterns  which  result  if  the  Aspect  Ratio  is 
increased  from  2  to  3.   It  is  shown  that  the  pattern  starts  to  split  in 
the  case  where  Re  =  10  and  the  Aspect  Ratio  =  3.   If  the  Reynolds'  number 
is  lowered  to  1,  the  split  becomes  complete.   This  latter  case  is  intar- 
esting  because  there  is  some  evidence  that  this  type  of  pattern  occurs 
in  Monterey  Bay.   It  had  been  speculated  that  the  splitting  of  the  cur- 
rent might  be  caused  by  the  presence  of  the  submarine  canyon.   But  as 
seen  later,  this  type  of  circulation  only  occurs  for  the  simple  cavity 
flow  model. 

All  of  the  vorticity  plots  are  similar.   A  typical  plot,  for  Reynolds' 
number  of  10  and  Aspect  Ratio  of  2,  is  illustrated  in  Figure  4.7.   It  can 
be  seen  that  the  walls  and  especially  the  corners  have  a  significant  effect 
on  the  vorticity.   It  is  doubtful  that  the  effect  of  the  sharp  corners  can 
be  easily  extrapolated  to  what  occurs  in  the  real  world. 

Defining  closed  circulation  as  that  which  occurs  in  any  area  where  the 
streamlines  form  a  closed  curve,  it  is  interesting  to  note  that  in  all 
cases  closed  circulation  does  occur.   It  will  be  seen  that  this  differs 
markedly  from  the  results  for  the  refined  model.   The  closed  circulation 
varies  from  occurring  in  only  the  corners  of  the  cavity  to  filling  the 
entire  cavity.   If  it  is  assumed  that  this  simplistic  model  approximates 
the  real  world,  then,  important  conclusions  about  a  bay  as  a  place  within 
which  to  deposit  effluent  can  be  deduced.   These  conclusions  are  discussed 
in  the  following  section. 
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D.   CONCLUSIONS 

As  was  seen  in  the  previous  section,  some  form  of  closed  circu- 
lation occurs  in  all  of  the  solutions  of  the  stream  function  obtained. 
If  this  is  the  case  in  reality,  then,  an  embayment  could  be  a  poor 
place  within  which  to  deposit  waste.   Should  the  effluent  be  deposited 
within  one  of  the  gyres,  diffusion  and  tidal  currents  would  have  to  be 
relied  upon  to  disperse  the  pollutant.   Diffusion  and  tidal  currents 
might  not  disperse  the  pollutant  at  a  rate  fast  enough  to  keep  the 
pollution  below  a  safe  level. 
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V.   REFINED  MODEL 

A.   GOVERNING  EQUATIONS 

The  development  of  the  final  equations  from  Equations  3.23  and  3.24 
lacks  only  the  formulation  for  the  bottom  friction.   The  bottom  friction 
terms  in  the  vorticity  transport  equation  (3.23)  are  not  in  a  form  which 
can  be  numerically  solved;  therefore,  seme  other  relationship  must  be 
found  from  which  they  can  be  computed.   Ihe  following  is  a  reasonable 
approach  to  the  problem  which  yields  a  form  for  the  bottom  friction  which 
can  be  solved  numerically. 

1.   Bottom  Friction  Equations 

This  development,  due  to  its  length,  is  derived  in  detail  in 
Appendix  A.   The  short  explanation  contained  here  will  introduce  the 
method  used  to  obtain  a  relationship  for  the  bottom  friction.   Consider- 
ing Ekman's  solution  for  slope  currents,  it  can  be  shown  that 

2 

d  w  _  ifw   g_  cU^ 

,2   A     A   9n 
dz     v     v 

where  the  complex  velocity  field  is  given  by 

w  =  u  +  iv 

The  solution   to   this   equation   is 

ig    /£Osh_a£        -  >     d^ 
f        cosh   ah  8n 


where 


a-   (r) 

V 


Now,  the  bottom  stress  equals 
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b   A  . 

T _  _V  QW  I 

ph  "  h   3z' 

z  =  -h. 


therefore 


£  =  ^  (f«)  tanh  ah  .£  (5.1) 

By  integrating  the  velocity  over  depth  it  follows  that  the  transport, 
Tr,  equals 

Tr  =  (|£)  (|^)  -  (tanh  ah  -  ah)  (5.2) 

i     dr\      a 

then  equating  5.1  and  5.2 

:L_  0    (11)    Tr    r      tanh   ah  i 

ph  h  tanh  ah   -  ah 

Reducing   this    equation   to   its   x,y   components   and  nondimensionalizing 

yields 

I —  =  Y    [-  —  V     +  —  V    ]  (5.3) 

pH  l      H       x        H        yJ  K        J 

xby  R?  R1 

I —  zzFr-^li}/     .illy]  (5#4) 

pH      l   H   x   H   yJ  v    J 

where  Rl  and  R2  are  coefficients  which  are  functions  of  ax  which  in 

turn  is  equal  to 

?  1/2 

/2  ,fh  \ 
ax  =  ^  (— ) 
v 

Assuming  that  the  depth  at  a  given  point,  once  specified,  remains  con- 
stant and  f  and  A  are  constant,  then  Rl/H  and  R2/H  need  only  be  computed 
once  at  each  grid  point  in  the  beginning  of  the  numerical  program  and 
stored  for  future  use  in  each  iteration  of  the  vorticity  equation.   A 
closer  look  at  the  equations  shows  that  the  bottom  friction  is  represented 
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by  depth-dependent  coefficients  multiplied  by  the  components  of  the 
velocity.   This  follows  directly  from  the  fact  that 

F  [-  Rl  r^  +   R2  t^-]    =   -   F  [Rl-V  +  R2'U] 
n        H 

The  other  bottom  stress  component  takes  a  similar  form. 

Table  II   contains  values  for  Rl  and  R2  which  are  computed  from 
the  full  equations.   The  frictional  coefficients  and  approximations 
used  to  calculate  them  for  large  ax  are  plotted  in  Figure  5.1  and  Figure 
5.2.   The  approximations  are  given  in  Appendix  A.   It  can  be  seen  from 
the  plots  that  the  approximations  are  excellent  as  long  as  ax  is  greater 
than  2.0.   Below  this  value,  the  full  equations  must  be  used.   Relating 
this  to  the  depth,  it  is  found  that  the  approximations  are  good  until 

h  <  3/A   (meters) 

At  very  small  values  for  the  argument  ax,  which  is  dealing  with  depths 
on  the  order  of  centimeters,  the  solution  oscillates  about  zero  and 
the  validity  of  the  solution  is  questionable.   This  is  not  investigated, 
because  at  depths  shallower  than  five  meters  the  surf  zone  is  encountered 
and  this  regime  is  controlled  by  a  different  set  of  processes. 
2.   Final  Form  for  the  Vorticity  and  Stream  Function 

Equations  for  the  vorticity  and  stream  function  of  the  refined 
model  were  derived  in  Section  III.   These  equations  with  slight  modifi- 
cations, depending  upon  which  forces  are  being  neglected,  are  used  in 
this  section.   Substituting  for  the  bottom  friction,  the  equations 
arrived  at  in  Section  III-l,  Equations  3.23  and  3.24  become 

V  Re"  V^  "  H  J(H/'Z)  "  H  J(VF'F)  "  (Z+2F)   J(HjVF) 

H  (5.5) 

U   H    x    H   V        H    x    H   V 

3    y  x 
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TABLE  II 
Frictional  Coefficients 


ax 

Rl 

R2 

.5 

-0.19995 

6.00305 

1.0 

-0.19899 

1.51133 

1.5 

-0.19509 

0.69132 

2.0 

-0.18587 

0.41530 

2.5 

-0.17090 

0.29412 

3.0 

-0.15284 

0.22951 

3.5   . 

-0.13526 

0.18872 

4.0 

-0.12014 

0.15995 

4.5 

-0.10773 

0.13849 

5.0 

-0.09756 

0.12197 

6.0 

-0.0819  7 

0.09836 

7.0 

-0.07059 

0.08235 

8.0 

-0.06195 

0.07080 

9.0 

-0.05517 

0.06207 

10.0 

-0.04972 

0.05525 

15.0 

-0.03325 

0.03563 

20.0 

-0.02497 

0.02628 

25.0 

-0.01998 

0.02082 

30.0 

-0.01666 

0.01723 

35.0 

-0.01428 

0.01470 

40.0 

-0.01250 

0.01282 

45.0 

-0.01111 

0.01136 

50.0 

-0.01000 

0.01020 

100.0 

-0.00500 

0.00505 

200.0 

-0.00250 

0.00251 

300.0 

-0.00167 

0.00167 

400.0 

-0.00125 

0.00125 

500.0 

-0.00100 

0.00100 

1000.0 

-0.00050 

0.00050 

2000.0 

-0.00025 

0.00025 
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Figure  5.1   Frictional  coefficient,  Rl 
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Figure  5.2   Frictional  coefficient,  R2 
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and 

Z  =  fj  V2Z  "  ^r  (Y  H   HH)  (5.6) 

H       R2   x  x    y  y 

Taking  the  indicated  partials  in  the  bottom  friction  term  gives 

ZT  =  Re  ^  "  H  J°',Z)  ~  H  J(4,>F)  "  ^f^  J(H^) 

H  (5.7) 

+  eF  [|^  V2^  +  J(§^,¥)  +  (I2-)  »  +  (f2-)   fx] 

y  x 

Tha  equations  for  the  refined  model  contain  the  volume  transport 
stream  function  rather  than  the  velocity  stream  function  which  is  con- 
sidered in  the  simple  cavity  flow  problem.   This  study  is  interested 
only  in  the  direction  of  the  flow  and  net  the  volume  transported.   The 
velocity  stream  function  has  the  property  that  the  direction  of  the  flow 
is  at  all  points  tangent  to  the  streamlines.   Here,  it  will  be  proven 
that  the  volume  transport  stream  function  has  the  same  property. 
Considering 

v-w  =  u|^+  v  I- 

dx     oy 
therefore 


V-VY  =  U  (+  UH)  +  V(-UH) 


thus 


v-vy  E  o 

This  means  that  the  velocity  is  perpendicular  to  the  gradient  of  the 
volume  transport  stream  function  which  in  turn  means  the  velocity  is 
everywhere  tangent  to  the  streamlines. 
3.   Convergence  and  Stability 

The  rate  of  convergence  for  the  refined  model,  as  expected,  is 
much  slower  than  for  the  simple  cavity  flow  problem.   To  obtain 
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convergence  to  within  one  percent  of  the  original  residual,  it  takes 
between  two  and  three  hours  of  computer  time.   On  the  average,  the 
equations  are  converged  until  the  residual  is  one-half  of  one  percent 
of  the  original  value.   The  proper  time  step  is  determined  experimentally 
and  is  found  to  be  about  one-fiftieth  of  the  time  step  used  in  the 
simple  cavity  flow  section.   Convergence  to  within  one  percent  of  the 
original  residual  required  approximately  600  time  steps. 

Although  the  Arakawa  Jacobian  is  a  stable  representation,  it  is 
found  that  the  bottom  topography  introduced  must  be  smoothed,  eliminating 
some  detail,  in  order  to  obtain  the  desired  degree  of  convergence. 

B.   INVESTIGATION  OF  THE  RELATIVE  IMPORTANCE  OF  THE  CORIOLIS  FORCE  AND 
BOTTOM  FRICTION 

One  question  which  occurred  during  the  development  of  this  study  was: 
how  important  are  the  Coriolis  force  and  the  bottom  friction?   To  inves- 
tigate this  question  the  nondimensional  parameter,  e,  and  the  non- 
dimensional.  Coriolis  force,  F,  is  examined  since  both  forces  contain 
these  terms. 

1.   Coriolis  Force 

The  magnitude  of  the  Coriolis  force  and,  to  a  degree,  the  bottom 
friction  depend  upon  the  magnitude  of  the  nondimensional  parameter,  e. 
Recalling  that 

e=f^  (5.8) 


where 


and 


2ft  cos  0 
3=  °- 


F  =  F  +  Y 
o 
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where 

Y  -  Direction  of  true  north 

F  -  t   tan  9 
o    L      o 

ft  ~  (7.29x10   rad/sec)  angular  rotation  of  earth 

o 

r  -  (6.38x10  cm)  radius  of  earth 

L  -  (1.609x10   cm)  characteristic  length  =  one  mesh  length 
in  x  direction 

0  =  (37°)  latitude  for  this  study 
o 

Carrying  out  the  indicated  arithmetic  operations  yields 

-13   -]     -1 

$  =  1.8  x  10    cm   sec 

and 

_   6.38  x  108  .    ,_o   _ 
F  =  c  tan  37  +  Y 

1.61  x  10 


which  reduces  to 


F  =  3  x  103  +  Y 


It  can  be  seen  from  5.8  that  once  the  characteristic  length  is  deter- 
mined, then  E  depends  only  upon  the  characteristic  velocity,  U  .   Fron 
the  definition  for  the  Reynolds'  number,  a  relationship  for  U   can  be 

determined. 

ReA 

u  =T-E 

Substituting  into  5.8  gives 

SL3 
e  = 


R6AH 

Assuming  a  value  of  10   for  A^  and  substituting  the  values  obtained  for 
3  and  L,  then  e  becomes 

7.64  x  10~5 
e  =  - 


Re 
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This  is  a  convenient  relationship  because  solutions  for  various 

Reynolds'  numbers  are  desired;  therefore,  once  the  Reynolds'  number 

is  specified,  the  value  for  e  is  determined. 

Assuming,  at  first,  that  the  Coriolis  force  is  constant,  the 

term  appealing  as  a  coefficient  in  front  of  both  the  bottom  friction 

and  topographic  tendency  terms  is 

eF  =  eF 
o 

Direct  substitution  yields 

.23 

eF  =  - — 
Re 

The  inclusion  in  the  model  of  the  constant  term,  eF,  is  necessary  and 

presents  no  problems.   Now,  the  effect  of  the  inclusion  of  the  variation 

of  the  Coriolis  term  with  Latitude  is  examined.   The  term  being  examined 

is  given  by 

|  J(Y,F)  =  |  [V  F  -fF] 
H  H   x  y    y  x 

It  has  been  determined  that 

eF  =  3  x  103e  +  Ye 
Thus,  if  it  is  assumed  that  the  grid  is  aligned  such  that  Y  is  in  the 
direction  of  True  North,  then 

(eF)     =  7.64  xlO"5 


but 


therefore 


y  Re 


(eF)      =   eF 

y         y 


§  [T  F  _  ^y0]  = 7-64  x  10~5  v* 


H    L   x  y         jf    x  Re  H 
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This   terra  is  at  least  several  orders  of  magnitude  less  than  the  other 
terras  in  the  equation  and  is,  therefore,  negligible.   Thus,  although 
the  variation  of  the  Coriolis  force  is  negligible  due  to  the  scale  used 
in  this  study,  the  presence  of  a  constanc  Coriolis  force  will  be  shown 
to  have  an  effect  upon  the  solutions  which  are  obtained. 
2.   Bottom  Friction 

Ex£mining  the  bottom  friction  terms  in  Equation  5.7  it  can  be 
seen  that  the  frictional  coefficients,  Rl  and  R2 ,  appear  in  all  of  the 
expressions.   Looking  at  Table  II  it  is  seen  that  the  frictional  coeffi- 
cients rapidly  become  small  and  the  change  of  the  coefficients  in  the 
horizontal  direction,  due  to  the  change  in  the  bottom  topography,  becomes 
negligible.   Considering  eF,  with  Reynolds'  numbers  greater  than  one, 
this  factor  which  is  multiplied  into  the  frictional  term  is  less  than 
.2.   Thus  at  Reynolds'  numbers  greater  than  one  or  at  depths  greater 
than  thirty  meters  the  bottom  friction  becomes  negligible.   Since  this 
model  uses  a  grid  spacing  of  0.5  miles  in  the  Y  direction  and  1.0  miles 
in  the  X  direction,  the  area  covered  and  the  topographic  changes  are 
large.   The  only  place  where  the  bottom  friction  would  have  any  signifi- 
cance is  at  the  few  grid  points  immediately  surrounding  the  boundaries 
where  conditions  are  only  approximated.   At  very  shallow  depths,  less 
than  five  meters ,  the  equations  no  longer  represent  the  current  regime 
because  the  surf  zone  is  encountered.   Also,  considering  that  the 
velocity  profile,  irrespective  of  bottom  friction,  is  generally  a  maximum 
at  the  surface  and  decreases  with  depth,  the  effect  of  bottom  friction 
is  minimized  in  this  situation  compared  to  the  situation  where  the 
velocity  is  constant  with  depth.   For  these  reasons,  it  is  concluded 
that  it  is  unrealistic  to  assume  that  the  inclusion  of  the  bottom 
friction  would  increase  the  accuracy  of  this  model. 
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C.   NUMERICAL  APPROACH 

1.  Method  for  Solving  the  Vorticity  Equation 

The  full  vorticity  equation,  Equation  5.7,  neglecting  the 

variation  of  the  Coriolis  force  with  Latitude  and  the  bottom  friction. 

reduces  to 

ZT  -  |^  V2Z  -  |  J(Y,Z)  -  (Z~^F)  J(H,¥)  (5.9) 

H 

in  finite  difference  form  this  becomes 

_n+l  _n   .    Zk    .+Zk+|  .-2Zk  .   Zk  .^+Zk+1     -2Zk  . 
Z =Z__  =  1_  r  i+l,j i-l ,  j i^    1,3+1   1,3-1   i,-| 

(Ax)  (Ay) 

(5.10) 

1  (Zi  i  +  eF) 

i,j  H2  . 

i,J 

where  the  Jacobians  in  the  above  equation  are  solved  by  Arak^wa's 
technique  which  is  outlined  in  Section  II-B-3.   For  programming  purposes, 
the  above  equation  is  manipulated  to  yield  a  similar  equation  which  mini- 
mizes number  of  computations  needed  to  obtain  a  solution. 

2.  Optimum  Overrelaxation  for  Solving  the  Stream  Function  Equation 
Substituting  the  appropriate  finite  difference  forms  into 

Equation  5.6  yields 

z  =  -1 t  1+1>J  i-i».i   LJ  +  1»-i+1  i».i-i  x»-i] 

Hi,j       (Ax)2  (Ay)2 

_    _i r         1+1,3         1-1  ».1  1+1,  .1  1-1,  .1  .  1,J+1         1,1-1  1,3+1         1,1-1      , 

2    L  2  2 

HT  .  4 (Ax)  4 (Ay) 

i  >  J 

2 
An  expression  for  ¥ .  .is  desired;  therefore,  first  multiply  by  (Ax) 

i » J 

and  then  collect  terms  obtaining: 
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2  2 

Av      k+1  Ay 

2  [1  +  (t21)  ]  V.  .  =  (T»    ,  +  ¥ ,  .  .)  +  (7^)   (V.  ...  +  T,  .  .) 
Ay     i,j     l+l, J     1-1,3     Ay     1,3+1    i,j-l 

-  (Ax)   H   Zn   -  (Ax)   [   1+1'J   1"1>-1    1+1»J   1"1»J 
i,i  i,j   H.  .  /  /.  \2 

1,3  4 (Ax) 

+   i»3+l   1,3-1     i>3+l   1,3-1  j 
4 (Ay)2 

AX 
Letting  B  =  —  and  employing  Optimum  Overrelaxation 

k+1    ,  .   k      RoDt     k      k+1      ?  k        k+1 

rv.   -  (1.0  -  Ropt)  v.  .  +   opu0  [H*    ..+¥.  .  .  +  BZC1\  ,._  +  Y.   .) 
i»3  i,3    2(1+B  )       ^J      ,J       1,J      3-»J_1 

/  t  k+l    X  , 

-  (Ax)2  H.  .Zn  .  -   i+l,3   i-l, J    i+l, J   1-1,3 

i,3  i,3  4H 

-1- ,  J 

E      k        k+1 

-  75 (*.  ..1  ~  *.  ,  -,)  (H.  ,,,  -  H.  .   )] 

4H   .    i,3+l    1,3-1    i,3+l    i,3-l 

J- ,  J 

Defining  the  constants  and  reorganizing  :'.n  order  to  minimize  computations, 
this  equation  is  used  to  calculate  the  stream  function. 
3.   Boundary  Conditions 
a.   Inflow 

In  this  model  the  stream  function  is  not  related  to  the 
velocity,  U,  as  in  the  simple  cavity  flow  but  to  the  volume  transport, 
UH;  therefore,  the  stream  function  is  a  function  of  both  the  depth  and 
the  velocity.   It  is  assumed  that  the  inflow  boundary  layer  is  completely 
dominated  by  the  bottom  profile  with  frictional  effects  being  negligible. 
This  being  the  case,  the  input  velocity  is  assumed  constant  and  the  stream 
function  boundary  layer  is  determined  directly  from  the  bottom  topography. 
Knowing  that 

Y  =  -UH 

y 
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then  one  procedure  for  determining  the  input  boundary  layer  is 

2Ay  l,j 

and 

vi  • ,1  =  ^i  •  t  -  2Ay  UH1  . 

This  procedure  for  determining  the  input  boundary  layer  is  found  to  be 
unsatisfactory  as  it  is  in  the  case  of  the  simple  cavity  flow.   To 
solve  the  problem  of  specifying  the  inflow  stream  function  boundary 
layer,  the  grid  is  extended  five  grid  spaces  upstream  and  four  spaces 
downstream  as  transitional  regions.   The  bottom  is  specified  as  linear 
at  the  artificial  inflow  boundary  gradually  changing  into  the  "true" 
bottom  profile.   When  this  is  done,  a  well  behaved  input  profile  can  be 
determined  by 


thus 


«F  =  -UH 

y 


yi 


Y  -  -  I     UH  dy 

y0 

but  U  is  constant  and  H  =  CY,  therefore 


finally 


yl 

Y   = 

-UC  J       Y   dy 

yo 

2 

yi 

«P   = 

-ucf- 

yo 

where  C  is   the  slope  of   the   artificial   inflow  bottom 

H  -  R\ 


C   = 


Y-Y1 
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It  follows  directly  from  Equation  3.11  that  the  vorticity 
is  zero  at  the  inflow  since  the  V  component  of  the  velocity  is  zero  and 
the  U  component  is  constant. 

b.   Lateral  Boundaries 

The  rectangular  shape  of  the  cavity  is  retained  in  the 
refined  model.   It  is  assumed  that  the  effect  of  the  bottom  topography 
will  simulate  the  effect  of  realistic  boundaries.   This  is  achieved  by 
using  the  actual  coded  bottom  but  specifying  a  minimum  depth  to  be  used 
whenever  the  depth  becomes  too  shallow  or  the  boundary  goes  overland. 
This  minimum  depth  is  used  until  the  wall  is  reached.   Diagramatically , 
this  is  illustrated  in  Figure  5.3.   Using  this  method,  the  problems  of 
programming  variable  boundaries  and  its  associated  large  increase  in 
computer  time  is  bypassed.   In  doing  this  the  boundary  condition  at  the 
wall  still  must  be  specified. 

The  stream  function  at  the  walls  is  still  equal  to  zero  since 
no  slip  conditions  are  applied  which  implies  that  the  velocity  at  the 
walls  is  zero.   The  vorticity  can  be  determined  in  the  following  manner. 
First,  consider  walls  alligned  in  the  X  direction.   At  the  wall,  V  =  0 
and 

wal1     *  wall  ^<all 
thus  o 

wall  ^h  ;  u2 

y     H 

and 

Z  11  =  h  * 
wall   H  yy 

Considering  Taylor's  expansion: 
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Figure  5.3  Grid  for  refined  model. 
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wall+1  yKall     y^\      in      2  yy   ,  /  y; 
^       -^   'wall  'wall 


therefore 

f   _   wall+1 
yy    (Ay)2 

Substitution  yields 

2T 
_  1_   wall+1 

wall  ""  H.  .  "   .  N2 
i,J   (Ay) 


Vertical  walls  are  handled  in  a  similar  manner  yielding 

2V 
1      wall+1 


wall   H.  .     ,.  N2 
i,j     (Ax) 

Sharp  corners  are  handled  in  the  same  way  they  are  handled 
in  the  simple  cavity  flow  problem.  It  is  expected  that  the  effects  of 
the  boundaries  on  the  solution  will  be  diminished  due  to  the  effect  of 
the  bottom  topography. 

c.  Downstream  Continuation 

The  same  method  as  the  one  employed  for  the  simple  cavity 

flow  problem  is  used.   This  method  gives 

w        =  y     =  2^       -  ¥ 
outflow    IL,J     IL-1,J     IL-2,J 

and 

outflow    IL,J    IL-1,J 

d.  "Lid" 

It  is  assumed  that  the  "Lid"  is  far  enough  from  the  cavity 
so  that  the  assumption  that  there  is  no  volume  transport  across  this 
imaginary  boundary  will  not  seriously  affect  the  flow  patterns  created. 
It  is  assumed  that  the  current  at  the  "Lid"  is  constant  and  parallel  to  the 
"Lid".   This  makes  the  "Lid"  a  stream  line  and  as  in  the  simple  cavity  flow 
problem: 
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Zi,LID  "  Zi,LID-l 
e.   Bottom  Topography 

Two  basic  bottom  topographies  are  used  in  this  study.   The 
first  bottom  topography  (Figure  5.4)  used  is  a  hypothetical  one  which 
is  very  simple  in  order  to  minimize  any  instability  it  might  generate. 
This  bottom  is  obtained  by  making  the  depth  at  the  "Lid"  a  constant  1600 
meters  and  assuming  a  linear  decrease  of  depth  along  the  inflow  boundary 
to  zero  at  the  first  grid  after  the  wall.   This  makes  the  value  at  the 
wall  the  minimum  non-zero  depth.   All  of  the  other  columns  exterior  to 
the  cavity  are  equivalenced  to  the  inflow  column.   The  entire  grid  is 
then  relaxed,  holding  the  boundaries  constant  until  a  steady  state  is 
reached.   The  presence  of  the  cavity  craated  the  bottom  topography  in 
Figure  5. 4. 

To  obtain  a  bottom  topography  which  simulates  that  of 
Monterey  Bay,  a  bathymetric  chart  of  the  bay  is  overlain  by  a  50  x  80 
grid  and  the  depths  are  coded  at  each  grid  intersection  and  then  placed 
on  computer  cards.   The  field  is  then  read  into  the  computer  and  stored 
for  use  in  the  program  as  needed.   The  maximum  depth  of  1600  meters  is 
chosen  because  it  is  approximately  the  maximum  level  of  no  motion  used 
in  geostrophic  current  calculations.   The  field  is  smoothed  by  averaging 
the  surrounding  four  grid  points  in  order  to  obtain  the  central  grid 

value ,  i.e.: 

H...  ,  +  H.  -  .  +  H.  .    +  H. 

H    =   i+l  >.1     r-1, .1     i,.i+l    1..1-1 

It  is  necessary  to  smooth  the  field  in  order  to  eliminate  inconsistencies 
in  the  bottom  topography.   In  order  to  retain  the  effect  of  the  canyon, 
the  boundaries  of  the  1600  meter  contour  are  preserved.   This  boundary 
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is  marked  by  the  dark  line  in  Figure  5.5.   Figures  5.6 — 5.9  show  a 
progression  from  the  original  field  to  one  smoothed  five  times.   It 
can  be  seen  that  after  only  five  smoothing  passes  only  the  general 
outline  of  the  canyon  is  still  preserved. 
4 .   Computational  Sequence 

The  computational  sequence  for  this  portion  of  the  study  is  very 
similar  to  the  one  used  in  the  simple  cavity  flow  section  with  a  few 
additions.   The  first  difference  is  that  the  bottom  topography  must 
either  be  generated  within  the  program  or  must  be  read  into  the  program 
from  some  external  source:   tape,  disc,  cards,  etc.   The  other  main 
difference  is  that  the  vorticity  transport  equation  is  much  more  complex 
and  some  of  the  dimensioned  variables  occur  more  than  once  in  the 
equation.   In  order  to  minimize  computer  time,  at  each  time  step,  the.se 
multiple-subscripted  variables  are  equivalenced  to  non-subscripted 
variables  before  the  vorticity  transport  equation  is  solved. 

D.   RESULTS 

The  vorticity  equation  used  for  the  refined  model  is  given  by 

ZT  =  ^  V2Z  -  I  J(^Z)  -  (Z+e2Y)    J(H,Y)  (5.11) 

H 

In  the  above  equation  the  stream  function  is  the  volume  transport  stream 

function.   There  are  two  cases  which  develop  naturally  from  this  equation, 

Case  I  occurs  when  eF  equals  zero.   Physically,  this  is  the  situation 

that  arises  when  the  Coriolis  force  is  neglected.   When  eF  is  assumed  to 

be  a  constant  other  than  zero,  then  only  the  change  of  the  Coriolis  force 

with  latitude  is  being  neglected.   This  is  case  II.   Solutions  for  both 

cases  are  obtained  and  discussed  in  the  next  two  sections.   A  measure  of 
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the  importance  of  the  various  terms  in  Ecuation  5.11  can  be  obtained  by 
varying  the  Reynolds'  number  and  eF. 
1.   Neglecting  the  Coriolis  Force 

For  this  case,  results  are  obtained  for  both  bottom  topographies 
discussed  in  Section  V-3—e.   In  the  solutions  obtained  here  the  flow  is 
from  the  left  to  the  right  representing  flow  from  the  South  to  the  North. 
For  the  hypothetical  bottom  topography  illustrated  in  Figure  5.3,  due  to 
its  symmetry,  the  solutions  are  the  same  for  flow  in  either  direction. 
This  is  true  only  because  the  input  velocity  is  assumed  constant  and 
therefore  the  input  vorticity  is  zero,   "his  is  tested  in  several  computer 
runs  by  reversing  the  position  of  the  inflow  and  outflow  and  differencing 
in  the  opposite  direction.   No  differences  in  the  solutions  develop  from 
reversing  the  direction  of  the  flow,  but  it  is  expected  that  if  gyres 
occur  which  are  not  in  the  center  of  the  cavity  then  the  two  solutions 
would  have  the  gyre  positioned  on  opposite  sides  of  the  cavity.   If  there 
is  an  input  velocity  boundary  layer,  different  solutions  would  be  obtained 
for  flows  in  opposite  directions  because  of  the  change  in  the  sign  of  the 
input  vorticity. 

Using  the  hypothetical  bottom  topography  illustrated  in  Figure 
5.3,  the  volume  transport  stream  function  plot  in  Figure  5.10  is  the 
result  obtained.   Computer  runs  for  Reynolds'  numbers  of  .01,  .1,  1,  10, 
100  are  made  and  all  of  the  results  are  qualitatively  identical  to  this 
plot.   The  plot  of  the  volume  transport  stream  function  obtained  when 
the  bottom  topography  simulating  Monterey  Bay  is  used  is  contained  in 
Figure  5.11.   The  results  for  this  bottom  topography  are  also  identical 
over  the  range  of  Reynolds'  numbers  from  .01  to  100. 

It  can  be  seen  from  the  results  for  the  hypothetical  bottom 
topography  that  within  the  bay,  the  streamlines  follow  the  bottom  contours. 
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This  means  that  the  direction  of  the  flow  within  the  bay  is  generally 
along  the  bottom  contours;  but  at  the  larger  depths  encountered  outside 
of  the  bay,  the  direction  of  the  flow  does  not  follow  the  contours.   The 
fact  that  trie  flow  does  not  follow  the  contours  outside  of  the  bay  while 
it  does  inside  the  bay  is  undoubtedly,  not  only,  the  result  of  shallower 
depths  within  the  bay,  but  also,  due  to  the  fact  that  the  volume  trans- 
port in  the  bay  is  much  slower  than  outside  of  the  bay.   For  the  bottom 
topography  simulating  Monterey  Bay  the  entire  cavity  is  filled  with  a 
gyre.   Again,  no  significant  difference  is  discerned  between  the  solutions 
obtained  over  the  range  of  values  investigated.   It  is  also  noticed  that 
the  cavity  seems  to  cause  a  disturbance  downstream  from  the  cavity. 
Although  the  field  downstream  is  not  long  enough  to  be  absolutely  certain, 
it  is  reasonable  to  assume  that  this  is  a  wave  induced  into  the 
flow  by  the  bottom  topography.   For  comparison  purposes,  the  simple 
cavity  flow  problem  is  solved  for  the  case  of  Reynolds'  number  equal  to 
.01  and  Aspect  Ratio  of  2.   The  result  is  contained  in  Figure  5.12.   It 
is  evident  what  a  significant  difference  the  inclusion  of  the  bottom 
topography  has  upon  the  resultant  solution. 
2.   Assuming  a  Constant  Coriolis  Force 

After  examining  the  case  where  the  Coriolis  force  is  neglected, 
the  effect  of  including  a  constant  Coriolis  force  is  investigated.   A 
range  of  values  of  the  Reynolds'  number  and  Coriolis  force  must  be 
investigated  because  these  parameters  depend  upon  the  characteristic 
scales  used.   Since  some  of  the  scales,  such  as  the  length  scale,  cannot 
be  determined  exactly,  solutions  must  be  obtained  for  a  range  of  values 
of  the  nondimensional  parameters.   It  is  then  assumed  that  the  solution, 
which  models  the  true  circulation  the  best ,  falls  within  this  range  of 

values. 
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Using  the  hypothetical  bottom  topography  shown  in  Figure  5.3, 
solutions  are  obtained  for  eF  equal  to  .23  and  23.   These  are  contained 
in  Figures  5.13  and  5.14,  respectively.   Each  figure  represents  solutions 
for  Reynolds'  numbers  of  .01  and  1  since  the  solutions  do  not  differ  over 
this  range  of  values.   As  can  be  seen  a  change  of  two  orders  of  magnitude 
of  eF  has  some  effect  on  the  flow  pattern  generated;  whereas,  a  change 
of  two  orders  of  magnitude  of  the  Reynolds'  number  has  a  negligible 
effect  upon  the  solution.   This  means  that  over  this  range  of  values  the 
effect  of  the  Coriolis  force  interacting  with  the  bottom  topography, 

^  j(h,<o 

H 
has  a  larger  effect  upon  the  solution  than  the  frictional  term, 

F"2z 

Re 
Still,  no  closed  circulation  occurs  in  any  of  the  solutions  obtained  for 
the  hypothetical  bottom  topography. 

The  next  step  is  to  obtain  solutions  for  various  values  of  the 
Reynolds'  number  and  the  Coriolis  force  for  the  bottom  topography  simu- 
lating Monterey  Bay,  Figure  5.4.   Solutions  are  obtained  for  the  follow- 
ing values  of  the  parameters : 

eF  Re  Figure 

.23  x  10"2  .01,  1  5.15 

.23  x  102  .01,  1  5.16 

1.00  x  103  .01,  1  5.17 

It  is  immediately  obvious,  that  with  the  Monterey  Bay  bottom  topography, 
closed  circulation  occurs  for  all  of  the  solutions  obtained.   Also,  the 
solution  for  the  hypothetical  bottom  fails  to  yield  closed  circulation 
over  the  same  range  of  values  for  the  two  nondimensional  parameters,  Re 
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and  eF.   It  can  be  concluded,  then,  that  the  closed  circulation  is  the 
result  of  the  presence  of  the  Monterey  Submarine  Canyon  and  that  the 
bottom  topography  has  the  most  significant  effect  upon  the  circulaticn 
pattern  generated. 

E.   CONCLUSIONS 

The  foundation  for  the  development  of  a  good  model  for  investigating 
the  current-driven  circulation  in  an  embayment  has  been  established. 
Included  in  this  model  are  the  effects  of  advection,  planetary  vorticity 
tendency,  topographic  vorticity  tendency  and  lateral  shear  stress.   The 
relative  order  of  magnitude  of  the  bottom  friction  and  the  change  of 
Coriolis  force  with  latitude  is  investigated  and  determined  to  be 
negligible:  compared  to  the  other  acting  forces  and  assumptions  that  are 
made  in  the  development  of  this  first  stage  of  the  model. 

There  are  two  areas  in  particular  where  the  model  could  be  improved. 
First,  effects  of  the  inclusion  of  the  bottom  friction  should  be  investi- 
gated further.   Effects  of  the  bottom  friction  diminishes  rapidly  with 
increasing  depth,  but  in  shallow  areas  around  the  boundaries  and  within 
the  bay,  the  bottom  friction  might  have  a  significant  influence  upon  the 
circulation  patterns.   Second,  some  of  the  boundary  conditions  should  be 
refined.   The  walls  representing  the  coastline  might  be  represented  in 
a  more  realistic  manner.   This  would  go  hand  in  hand  with  the  inclusion 
of  the  bottom  friction.   Also,  the  inflow  and  outflow  should  be  represented 
by  less  restrictive  equations. 

The  representation  for  the  bottom  friction  is  an  approach  suggested  by 
Dr.  Jerry  Gait  and  developed  in  this  thesis  by  the  author.  It  is  believed 
to  be  a  new  approach  and  one  worth  investigating  further.   This  approach 
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is  of  special  interest  because  the  bottom  friction  is  a  linear  function  of 
the  velocity  making  it  amenable  to  analytic  solutions. 

It  is  assumed  in  this  model  that  the  input  vorticity  is  zero;  because 
of  this,  the  solution  does  not  depend  upon  the  sign  of  the  vorticity  gen- 
erated at  the  boundary  which  is  a  function  of  the  direction  of  the  flow. 
If  a  velocity  boundary  layer  exists,  as  many  studies  propose,  then  the 
vorticity  will  not  be  zero  at  the  inflow  and  the  solution  will,  in  addition 
to  the  other  variables,  become  a  function  of  the  direction  of  the  flow. 
This  is  a  case  which  should  be  investigated  not  only  to  find  out  the  effect 
of  a  boundary  layer  upon  the  solution  but  to  determine  the  effect  on  the 
solution  of  changing  the  direction  of  the  flow. 

Results  of  this  study  indicate  that  the  occurrence  of  closed  circu- 
lation in  Monterey  Bay  is  a  distinct  possibility.   This  investigation 
indicates  that  the  presence  of  the  submarine  canyon  seems  to  increase  the 
possibility  of  the  occurrence  of  closed  circulation. 

In  conclusion,  it  can  be  stated  that  this  study  indicates  that  the 
bottom  topography  is  the  controlling  factor  in  determining  what  type  of 
circulation  pattern  is  to  be  expected  in  a  given  bay.   Over  the  ranges  of 
values  examined  for  Re  and  cF,  the  change  of  the  vorticity  from  one  case 
to  the  next  had  little  effect  upon  the  solution  of  the  stream  function. 
This  is  because  in  deeper  water  the  volume  transport  becomes  much  larger 
than  the  vorticity  and,  therefore,  the  stream  function  equation  becomes 

I  V2Y  _  I_^  (H/  H  +  y  H  )  =  o 
H       R2   x  x    y  y 

which  is,  for  all  practical  purposes,  independent  of  the  vorticity.,  This 
puts  additional  constraints  upon  the  types  of  flow  patterns  which  may  occur, 
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APPENDIX  A 

DEVELOPMENT  OF  BOTTOM  FRICTION 
Assuming  that  the  coefficient  of  vertical  eddy  diffusivity,  A  ,  is 
constant,  the  bottom  stress  is  proportional  to  the  change  of  the  mean 
velocity  gradient  with  depth.   Using  Ekman's  solution  for  slope  currents 
an  expression  for  the  velocity  as  a  function  of  depth  can  be  obtained,  but 
in  this  relationship  the  velocity  is  also  a  function  of  the  slope  of  the 
sea  surface.   In  order  to  eliminate  the  sea  surface  slope,  which  is  not 
easily  obtainable,  the  velocity  can  be  integrated  over  depth  to  give  an 
expression  for  the  volume  transport.   Equating  these  two  relationships 
an  expression  for  the  bottom  stress  which  does  not  depend  upon  the  slope 
of  the  sea  surface  can  be  obtained.   The  first  step  in  this  development 
is  to  obtain  the  analytic  solution  for  the  velocity  from  the  following 

three  equations 

2 
fv  +  A  1-a  .  -  1 1£  (A-!) 

V  3Z2     p  dx 

2 

-fu  +  A  ^--if  (A-2) 

v  8z2     p  9y 

P  =  -pgh  (hydrostatic  pressure)  (A-3) 

These  equations  can  be  derived  from  equations  II-l  and  II-2  with  the 
addition  of  the  following  assumptions 

1.  Steady  state  conditions  prevail 

2.  Non-linear  terms  are  negligible 

3.  Lateral  turbulent  stress  is  negligible 

Ekman  solved  these  equations  analytically  in  the  following  manner: 
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and 


3P  3h 

"3x~  =  "Pg  al 


3P  ah 

a?  =  ~pg  ^ 


Substitution  yields 


2 

8   u  _       Jfv  £_  _3h 

-    2            A  A     8x 

3Z                         V  V 


2 
9   v  _  j[u       j*_  ah 


thus 


defining 


it   follows 


.  .2        A  A      ay 

3z  v  v     J 


—2    (u+xv)    =  -  (u+iv)   +  —  (—  +  i  — ) 
3z  v  v  y 


w  =   u  +  iv 


3ri        3x  3y 


A_w  =  ±  fir +i_l£  (A_4) 

,2        X  A         A     3n  kA  q; 

dz  v  v 


Ekman   found   that   the   solution   to   this    equation  is 


w  i£       cos£_az  _  |c 

f  cosh   ah  3n 

where 

a=   (r) 

V 

Substituting  A-5  into  equation  A-A  readily  shows  that  it  is  a  solution. 

Given  an  expression  for  the  velocity,  expressions  for  the  bottom 
stress  and  the  transport  can  now  be  found.   The  bottom  stress  is  given  by 
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since 


b   A   - 

T _V  3w  I 

ph  =  h   3z  ' 

z=-h 


3w 


3? 


|5  =  a  (-ifi)  tanh  ah  , 
3  z      f  3n 


therefore 


b   A  a  . 

—  =  -t—   (rf)    tanh  ah  — - 
ph    h    f  3n 


(A-6) 


Determining  the  transport  is  as  easy  as  the  determination  of  the  bottom 


friction.   Knowing  that 

o 


then 


thus 


T  =  J  wdz 
-h 


•T   -  f  (i&)  (f )  (£2^i  -  1)  dZ 
r       J       f  3n   cosh  ah 


T  -  (^)  (|5")  [ ^— r  -1-  sinh  az 

r    f    3n   cosh  ah  a 


-  Z 


-h 


-h 


finally 


T  =  (§*)  (|£)  (-)  (tanh  ah  -  ah) 
r    r     dri   a 


Equating  A-6  and  A- 7  gives 


—  a  T  =  7—  a  (rf2-)    (— )  tanh  ah  [1  -  - — - — - 
h      r   h     f    3n  tanh  ah 


(A-7) 


from  this  the  desired  expression  for  the  bottom  friction  is  obtained 


b    A   n 
T     v  2  _   , 


tanh  ah 


] 


(A-8) 


This  form  is  not  convenient  for  programming.   Transforming  equation  A-8 
into  a  more  convenient  form  involves  separating  it  into  its  x,y  components 
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This  is  accomplished  in  the  following  manner:   first,  the  transport  is 
defined  in  terms  of  the  stream  function 


but 


T  =  T  +  i  T 
r    x      y 


T  =  UH  =  -Y 
x  y 


T  =  VH  =   ¥ 

y       x 


thus 

T  =  -¥  +  i  ¥ 
r    y      x 

secondly,  the  hyperbolic  terms  must  be  separated  into  its  components. 

Defining 

9    1/2 

,    .1/2  ,fhA       .1/2 
ah  =  i    (— — )    =i   x 

v 
therefore 


1/2 


T  ..  if  „   r    TANK  i  '  x 


where 


ph  '  h  Tr  *„  .1/2    .1/2  ]  (A  9) 

M  TANH  l    x  -  x    x 


.1/2        ._o    .   .   ,_o    ft   .    .    Jl 

l    =  cos  45+i  sin  45   =  — r-  +  i  — r  =  a+ia 


Substituting  for  the  arguments  of  the  hyperbolic  terms  yields 

1/2 

TANH  i  '    x  TANH  (ax+iax) 


_AMU  .1/2     .1/2    TANH  (ax+iax)  -  (ax+iax) 
TANH  l    x  -  l    x 

through  fundamental  hyperbolic  trigonometric  identities  and  some  simple 

complex  variable  manipulations  the  following  relationship  can  be  arrived 

at 

TANH  (ax+iax) =  Rl  +  i  R2 

TANH (ax  +iax)  -  (ax+iax) 

in  which : 
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R1      G-C 

G  +  2  [A-C] 


R2  = 3 

Z   G  +  2  [A-C] 
and 

A  =  ax  [COSH  2ax  +  COS  2ax] 
C  =  SINH  2ax  +  SIN  2ax 
D  =  SINH  2ax  -  SIN  2ax 


_    (COSH  2ax  -  COS  2ax) 
G  — 


(A-10) 


ax 


substituting  into  equation  A-9  gives 


b 

V  =  7T-   (-V   +iY  )  (Rl  +  i  R2) 

ph   h     y   x 

expanding  terms  yields 

bx      by    f 

■V-  +  i  "tt-  =  r  [-R1Y  +  R2^  +  i  (-R2Y  -  R1Y  )  ] 
ph      ph    h      x      y  x     y 

equating  real  and  imaginary  terms  leads  to  the  desired  bottom  friction 

components,  i.e.: 

bx 

^7—  =  -  (-R1Y  +  R2^F  ) 
ph    h      x      y 

and 

T^     f 

-^r—  =  -  (-R2T   -  RTF  ) 
ph    h      x      y 

Approximations  for  the  frictional  coefficients,  Rl  and  R2,  can  be  derived 

when  2ax  is  greater  than  4  by  eliminating  negligible  terms  in  the 

equations  A-10.   The  approximations  can  be  shown  to  be 

1  -  ax 


Rl  = 


2 
2 (ax)   -  2ax  +  1 
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R2  =  f 

2  (ax)   -  2ctx  +  1 

After  these  equations  were  developed  it  was  determined  from  the 

relative  order  of  magnitude  of  the  terms  involved  that  the  bottom 

friction  is  negligible;  therefore,  this  representation  was  not  tested. 

Plots  of  the  frictional  coefficients  and  tables  of  it's  values  are 

contained  in  Section  V.   By  way  of  comparison,  though,  it  can  be  shown 

that  the  results  from  this  method  can  be  interpreted  to  agree  with  the  values 

for  the  bottom  friction  used  by  Hansen  who  represented  the  bottom  friction 

by 


b    t         I   2a.  2      /TT^ 
t   =  (ru  vu  +v  ,  rv  vu  +v  ) 

where 

r  =  3  x  10~3 

This  representation  is  the  same  as 

b   0   ,„-3   .,2 
t  =  3  x  10   x  V 

where  V  is  the  velocity  vector.   The  method  described  in  this  thesis  gives 

approximately  the  same  value  for  "r"  over  the  depth  range  of  2000  -  3000 

meters  if  it  is  assumed  that  A  =  100.   A  big  difference  is  that  the 

v 

2 
coefficient  in  this  study  is  multiplied  by  V  instead  of  V  ,  but  if  it  is 

assumed  that  the  mean  currents  over  depth  investigated  are  of  the  order  of 

magnitude  of  1  cm/sec  then  the  results  from  both  methods  are  similar.   It 

is  impossible  to  compare  the  two  exactly  because  the  methods  for  computing 

the  bottom  friction  are  dissimilar.   This  method  does  allow  the  bottom 

friction  to  be  represented  in  a  rational  manner  and  has  the  added 

attraction  of  being  linear  in  velocity  making  it  amenable  to  inclusion 

in  analytical  solutions. 


Walter  Hansen,  Hydrodynamical  Methods  Applied  to  Oceanographic  Problems, 
Proceedings  of  the  Symposium  on  Mathematical-Hydrodynamical  Methods  of 
Physical  Oceanography,  p.  25-34.   Institut  fur  Meereskunde  der  Universitat 
Hamburg  (1962) . 
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APPENDIX  B 


SOLUTIONS  FOR  THE  SIMPLE  CAVITY  FLOW  PROBLEM 
This  appendix  contains  Calcomp  plots  of  solutions  to  the  simple  cavity 
flow  problem.   There  are  twenty-four  plots,  twelve  of  the  stream  function 
and  twelve  of  the  vorticity.   Two  different  techniques  are  used  to  obtain 
these  solutions.   Upwind  differencing  is  used  to  obtain  as  many  solutions 
as  possible  but  when  this  technique  failed  to  converge  the  second  method, 
employing  the  Arakawa  Jacobian  for  solving  the  advective  terms,  is  used. 
The  following  is  a  breakdown  of  the  technique  used  to  obtain  each 
solution. 

Reynolds '  Number 
1,10,100,1000 
1,10,100 
1000 
1,10,100,1000 


Aspect  Ratio 
1 
2 
2 
3 


Technique 

Upwind  differencing 

Upwind  differencing 

Arakawa 

Arakawa 
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